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The detection of gravitational waves from coalescing compact binaries would be a computationally 
intensive process if a single bank of template wave forms (one step search ) is used. In an earlier 
paper we had presented a detection strategy, called a two step search, that utilizes a hierarchy of 
template banks. It was shown that in the simple case of a family of Newtonian signals, an on-line 
two step search was ~ 8 times faster than an on-line one step search (for initial LIGO). In this 
paper we extend the two step search to the more realistic case of zero spin post 1,5 -Newtonian wave 
forms. We also present formulas for detection and false alarm probabilities which take statistical 
correlations into account. We find that for the case of a post 1,5 -Newtonian family of templates and 
J — , signals, an on-line two step search requires ~ 1/21 the computing power that would be required for 

■ the corresponding on-line one step search. This reduction is achieved when signals having strength 

G\ ' S = 10.34 are required to be detected with a probability of 0.95, at an average of one false event 

per year, and the noise power spectral density used is that of advanced LIGO. For initial LIGO, 
the reduction achieved in computing power is ~ 1/27 for S = 9.98 and the same probabilities for 
detection and false alarm as above. 
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I. INTRODUCTION 



A radiation reaction driven inspiral of a binary composed of compact massive objects (Neutron Stars, Black Holes) 
would emit gravitational waves that would lie, during the last few minutes before merger, in the sensitive bandwidth 
of Laser interferometric detectors like LIGO 0, VIRGO and GEO600. Even though most such events will produce 
a signal amplitude well below the noise rms at any given instant, the predictability of such wave forms would allow 
^C) • the use of pattern matching techniques like matched filtering to considerably improve their chances of detection ||. 

In matched filtering, the detector output is passed through a filter that is matched to the expected signal wave form 
in some optimal sense. If the maximum of the output crosses a pre-determined threshold, a signal is declared to be 
present in the data with a time of arrival given by the location of the maximum. The filtering of the detector output 
can, of course, be substituted with a cross-correlation against a suitable template wave form that is matched to the 
expected signal. 

Generically, the wave form from an inspiraling binary is an amplitude and phase modulated sinusoid both of whose 
instantaneous frequency and amplitude increase as the two bodies proceed towards merger. The signal becomes 
^ . "visible" in the output of a detector when its instantaneous frequency exceeds the lower frequency cutoff of the 
output's bandwidth. This moment can be taken as the time of arrival of the signal at the detector. Such a cutoff is 
required, for instance, in the case of ground based detectors because of excessive seismic noise at low frequencies. There 
would, of course, be an unknown phase offset at the time of arrival. In addition, the signal would be characterized by 
the masses and spins of the two components, the distance to the binary and geometrical factors like the direction to 
the binary and the orientation of the orbital plane. 

Thus, if matched filtering is used for such signals, it would be necessary to employ a bank of filters (or template 
wave forms) corresponding to different values of the signal parameters mentioned above. One would then compare the 
maximum over all the filtered outputs with a threshold. This strategy is usually called a one-step search. Even though 
the time of arrival, initial phase and distance can be handled easily, a one-step search would still be a computationally 
expensive proposition. Estimates |^],|5| of the computational power that is required for an on-line one-step search 
using post-Newtonian templates turn out to be ~ 200 Gflops (Giga flops) or higher. This is with the omission of 
various other signal processing overheads which can be expected to be present in a realistic situation. Therefore, it is 
desirable to have computationally less expensive detection strategies without, however, compromising too much on the 
performance afforded by matched filtering. One such strategy, called a two step hierarchical search, was investigated 
by us in an earlier work We found that this strategy reduces the computational cost of detection significantly 
without losing out on the performance of a one step search. This is because a two-step search utilizes information that 
was present in a one-step search but which was neglected. Namely, the correlation between templates which allows a 
coarse scan of the parameter space to predict the location of a threshold crossing peak among the filtered outputs. 

Our analysis was restricted in |6| (henceforth referred to as MD96) to the case of Newtonian wave forms and the 
noise power spectral density used was that of the initial LIGO. The main result of our analysis was that a two step 



1 



hierarchical search is ~ 8 times faster than the corresponding one step search. This gain was achieved when the 
detection probability desired was 0.95 for a signal to noise ratio of 8.8<r at an average rate of false events of 1/y. The 
Newtonian template family will, however, not be good enough for the detection of the true signal wave form 0. It 
was chosen in MD96 in order to keep the analysis simple since that work was in the nature of a first estimate and 
several other issues needed to be highlighted. This paper is an extension of the two step search to a more realistic 
family of signals and templates, namely, the post 15 -Newtonian family. 

Our choice is motivated by the result of Apostolatos || that a post 1 ^-Newtonian template family, having spin 
parameters (3 — a = 0, is adequate for the detection of signals up to (and possibly beyond) post 2 -Newtonian order, even 
for maximally spinning systems. However, this holds for spins that are aligned with the orbital angular momentum. 
In general, the signal from a misaligned system would suffer significant phase and amplitude modulations that can 
considerably reduce the SNR, for some detector-binary geometries, if non-spinning templates are used. The larger 
the opening angle between the orbital and total angular momentum, the larger is the fraction of detector-binary 
geometries which would be lost. But for moderate opening angles (~ 25° or less), a sizable fraction of signals can still 
be detected with the non-spinning templates. Therefore, a non-spinning post ^-Newtonian template family appears 
to be a realistic one to use. We choose our family of signals also to be the same since, as mentioned above, even higher 
order signal wave forms may be detectable using this family of templates. This choice of templates and signals should 
provide a realistic model for the assessment of a two step hierarchical search while keeping the analysis relatively 
simple. In the following we will refer to non-spinning post 15 -Newtonian wave forms as simply post ^-Newtonian 
ones. 

The main result of this paper is summarized in Tables | and Jolln Table |, the noise power spectral density (p.s.d.) 
expected for the initial LIGO |8| has been used while in Table |ll| the p.s.d. used is that which is expected for the 
advanced LIGO |J. Columns seven and eight of each table show the computational power required for an on-line 
two-step search (C^] inc ) and an on-line one-step (Celine) respectively for the same performance parameters. That 
is, a detection probability of 0.95 for all signals having a strength S lri \ n (as given in the captions of the tables) and 
an average false alarm rate of 1 false event/y. These values of the computational requirement have been obtained for 
various lengths of the input data segment which are tabulated in the first column. £ max is the length of the template 
having the lowest values for the binary masses, m\ and 7712, which we choose to be mi = 7712 = 0.5 M Q . The highest 
masses that we have used in our analysis are mi = to 2 = 30.0 Mq. These results show that a two step hierarchical 
search can reduce computational requirements by about a factor of ~ 25 in a realistic scenario. 

In MD96, the formulae used for detection and false alarm probabilities used the assumption of statistical inde- 
pendence between certain random variables. This assumption fails when templates are placed very closely and we 
were, thus, limited from exploring the small spacing case more thoroughly. In the present paper, we present a much 
improved formula for the detection probability that reproduces the Monte Carlo results quite well. It also suggests an 
alternative approach to Monte Carlo simulations for parameter estimation which should be orders of magnitude faster 
than the conventional approach but further investigations in this direction are postponed to a later work. We also 
show, somewhat qualitatively, that the assumption of statistical independence is justified as far as the false alarm is 
concerned. Thus, the results that we have obtained can be considered to be fairly accurate within the approximations 
that have been made due to other reasons, like the non-trivial boundary of the space of interest and the location 
dependence of the intrinsic ambiguity function. 

The rest of the paper is organised as follows. In Section |n| we apply the method of maximum likelihood detection 
to the post 1 ^-Newtonian family of signals. This is the rigorous formalism behind the matched filtering algorithm 
mentioned above. We end this section with the derivation of a test statistic whose value determines the choice 
between detection and non-detection. In Section III, we study the probability distribution functions of the test 
statistic. Formulas for the detection and false alarm probability are derived. In Section IV the problem of optimaly 
placing templates in a one-step search is investigated. We obtain some approximations regarding the placement 
geometry, number of templates etc. Section [v| is devoted to the two-step hierarchical search and also contains the 
results of this paper. We conclude with Section |v|. 



II. MAXIMUM LIKELIHOOD DETECTION OF POST 1 5 -NEWTONIAN SIGNALS 

The method of maximum likelihood detection entails the maximization of the posterior probability p(x(t) | O) over 
the signal parameters O. Here, x(t) is a segment of the output of a detector and p(A \ B) denotes the conditional 
probability of A given B. Actually, one should maximize p(Q \ x(t)) because it is x{t) which is given but when our 
a priori knowledge of the frequencies with which various values of the parameters can occur is negligible, this is 
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equivalent to the maximization oip(x(t) | 0) . The maximum is then compared with a fixed threshold and a detection 
is announced if the threshold is crossed. Maximum likelihood detection also achieves, the highest average detection 
probability (averaged over 0) for a given false alarm probability. Actually, this statement is only true approximately 



but the approximation becomes better as the signal to noise ratio becomes larger |10 . 

If the detector noise is assumed to be a stationary Gaussian process, defined by a (one-sided) power spectral density 
S n (f), the above maximization reduces to the maximization, over 0, of the quantity 



,s n (f) yjJ w ' ' 2j_ oa S n (f) 

where a stands for the Fourier transform of the corresponding time domain function and s(t; 0) stands for a 
member of the signal family. This motivates the definition of an inner product, 

(u(t),v(t))= J^JL^u(f)v*(f), (1) 

and a corresponding norm, 

HI = [( u > u )f' 2 ■ (2) 

Thus, maximum likelihood detection involves the computation of, 



A = max 
e 



(x(t),s(t; ©))-i( s (i; e),s(t; 0)) 



(3) 



This quantity is known as the test statistic. The test statistic is then compared with a pre-determined threshold to 
decide whether the given x(t) contains a signal or not. It should be noted that the properties of stationarity and 
Gaussianity are only approximations for the noise that will be present in the interferometric detectors. However, these 
approximations are expected to be quite good and we will assume such a noise in the following. 

Some clarification should be made here regarding the meaning of a detector output in the context of interferometric 
gravitational wave detectors. The final output at the photodetector would contain the response of the detector (as 
calculated using the geodesic deviation equation) convolved with the detector's transfer function jll| (which depends 
on the way the detector is configured, i.e., the kind of recycling used, the reflectivities of mirrors etc.). This signal 
would be buried in noise that would be white Gaussian noise if photon shot noise were the only source of noise. So, 
strictly speaking, the detection strategy should be for the detection of this convolved signal in white noise. However, 
it is easy to prove that this is equivalent to the detection of the deconvolved signal in a noise which has a power 
spectral density S n (f) such that when this combination is "passed" through a noise free interferometer, the output is 
the convolved signal and white noise, as would be observed actually. Clearly, S n (f) should be inversely proportional 
to the modulus squared of the detector's transfer function. When the noise at the output is not white, as would be 
the case in practice, S n (f) would be the power spectral density of the actual noise divided by the modulus squared 
of the detector's transfer function. Henceforth, by the detector output we would mean the deconvolved output which 
would have noise with a power spectral density given by S n (f) and the signal would be just the bare response of the 
interferometer i.e., the relative strain produced in the two arms as computed using the geodesic deviation equation. 



A. The post 1 5 -Newtonian signal 

In the case of coalescing binary signals expressed up to the post^-Newtonian order (spin parameter a = 0), the 
signal parameters are, © = {A, t ai To, T1.5). The parameter A is the (nearly) constant part of the amplitude 
of the signal that takes into account the distance to the binary and the relative orientation of the detector and the 
TT coordinate frames. The rapid rise of power in the seismic noise towards lower frequencies would require that the 
detector output be band pass filtered with a cutoff at some low frequency (usually assumed to be 40 Hz for the initial 
LIGO and 10 Hz for the advanced LIGO). Similarly, a cutoff f c will also be required at the high frequency end because 
of a rise in photon shot noise. Usually, f c is taken as f c — 1000 Hz. Thus, loosely speaking, the signal wave form 
from an inspiraling binary would "start" in the output of the detector at the time when its instantaneous frequency 
crosses f a . This instant is called the time of arrival of the signal and is denoted by t a . The phase of the wave form 
at t = t a is denoted by The remaining parameters have the dimension of time and depend on the masses of the 
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binary components. They have been called as intrinsic parameters of the wave form in contrast to t a and $ which 
are known as extrinsic parameters. 

Actually, the post 15 -Newtonian wave form should be characterized by three intrinsic parameters, (tq, t\, T1.5). The 
subscripts of r denote the post-Newtonian order at which that parameter occurs. Thus, tq is the Newtonian chirp time 
characterizing the lowest order wave form obtained using the quadrupole formalism and the post 1 -Newtonian wave 
form would have tq and t\ as its intrinsic parameters. However, we have assumed the spins of the binary components 
to be zero and, hence, we are left with only two independent intrinsic parameters which we have chosen as ri.g and To 
because they can be easily inverted to obtain the masses. In the expression for the wave form, however, we retain n 
with the understanding that it is dependent on ri.g and To- In an approximate sense, tq + t± — ri.g can be taken to 
be the time left to the final merger of the binary, starting from t = t a . 

In the following we will deal with the restricted form of the post 15 -Newtonian signal in which post-Newtonian 
corrections are applied to only the phase of the signal. The restricted wave form at any post-Newtonian level is 
expected to be a good model for the correct wave form at the same level ||. The restricted post 15 -Newtonian signal 
is, 



h(t; 6) = Aa(t — t a ; t ) cos 



dt'f(t'; r ,r 1 . 5 ) + $ 



(4) 



where, 



a(t) = li- 



ra 



-1/4 



and f(t; To, ri.g), the instantaneous frequency of the signal, is given by an implicit equation, 



t = Tq + Ti - T1.5 - t x 



,-s/3 



T\X 



T1.5X 



-5/3 



(5) 



(6) 



where x = f(t; tq, Ti.g)// a . The chirp times are given by the following expressions (G = c 
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Tl ' 5 = 8^ 



M 



(7) 
(8) 
(9) 



where, M is the total mass of the binary, /i is the reduced mass, 77 = [ijM and A4 = (/^M 2 ) 1 / 5 is the chirp mass. 
We have chosen tq and ri.g as our independent parameters. In terms of these parameters, 
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(10) 

(11) 

(12) 



Note that if t\ were used instead of ri.g, then the inverse relations for ra\ and m,2 would be more complicated and 
would have to be solved numerically. 

In Fig. [l] and Fig. ^, we have shown the (mi, m 2 ) plane (the binary masses), 0.5M Q < mi,m 2 < 3OAf , mapped 
into the (ri.g, tq ) space for f a — 40 Hz (initial LIGO) and f a = 10 Hz (advanced LIGO) respectively. The boundaries 
of the region, which we call the space of interest following Apostolatos || , can be obtained easily. The equation for 
the curve AC, corresponding to mi = m,2 or M — 4/i, can be found by directly substituting for [i and M from (|ll| ) 
and Eq. dfj). To obtain the other two segments, AB and CB, we use the following expression, 



25 6 , , 8 / 3 

— TQ{ir]a) ' = 



(mi + m 2 ) 1/3 



TOl?Ti2 



(13) 
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For every value of To, the line of constant tq will intersect AB at a point where one of the masses, say m 2 , is 0.5M Q . 
Similary, for BC, one of the masses, say mi, would be 30.0M Q at the point of intersection. The point where AB and 
BC meet falls on the To = To (30.0 A/©, 0.5 M Q ) line. Thus, for values of To larger than this, one of the masses in the 
above equation can be set to 0.5M Q and the equation can be solved for the other mass to yield the value for T1.5 at 
the point of intersection. For smaller values of To, one of the masses should be set to 30M Q and T1.5 be obtained as 
before. Thus, given a value of To, the two limits of T1.5 can be computed. This allows the area of the space of interest 
A to be computed using a standard 2D Quadrature algorithm (we use D01DAF of the NAg library). The area of the 
space of interest thus obtained is, 

._( 50.174 sec 2 for initial LIGO , > 

\ 20389.542 sec 2 for advanced LIGO ' ^ ' 

We can write the R.H.S. of Eq. (g) as 

h(t; 0) = Ah (t-t a ; 7-0,7-1.5) cos$ + Ah w / 2 (t - t a ; 7-0,7-1.5) sin $ , (15) 

ho{t; 7-0,7-1.5) = a(t; t )cos (J dt' f(t'; t ,ti. 5 )^ , (16) 

K/2{t;TQ,n.s)=a(t;To)cos(j dt' f(t'-, T ,n. 5 ) + ~\ . (17) 

This representation will be useful in what follows. 

The Fourier transform of h(t; 9) can be computed using the stationary phase approximation. For our purpose, it 
suffices to give its overall form here, the details being available in other sources [BUM]. For / > 0, 



h(f; 9) ex AT 7/6 exp[-2 7 r l (/< a + £ n^f)) + i$] , (18) 

i 

where the index i £ {0, 1, 1.5}. The functions tpi are independent of the signal parameters. For / < 0, the transform 
is constructed using the Hermitian property of the Fourier transform, h(f) — h*(—f), since hit, 9) is a real function. 



B. The test statistic and its computation 

Following the brief outline given earlier, the maximum likelihood detection strategy for post^-Newtonian signals 
would consist of the computation of a test statistic A given by, 

A = max [( x(t), h{t; 9)) - l/2( h(t; 9), h(t; 9) )] . (19) 
For the sake of convenience in the following, we adopt the notation 

0' = (*a,Tl.5,7-o), (20) 

= (n.5,7- o ). (21) 

Occasionally, we will also break up 9' as t a and 6. The maximization over the parameters A and $ can be carried 
out analytically to yield, 

. _ {x, qo) 2 + (x, q^/2) 2 - (x, q )(x, q 7r / 2 )(qo, 9^/2 ) , 00 , 
— max , 

c 1 - {qo, qir/2) 2 



where, 



q (t; 9) =N - 1 h Q (t- 9) , (23) 

q,/2(t; 9) =M-J- 2 K /2 {t; 9) , (24) 

Af - \\h {t; 9)\\ , (25) 

K/2 = \\K /2 (t; 9)\\ . (26) 
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No and ■Mr/2 are normalization constants that are chosen as above for later convenience. Since t a occurs in the phase 
of the Fourier transforms of ho(t — t a ; 9) and h^/ 2 (t — ta', 9), No and N„n are independent of t a . We call qo(t;0) 
and q^/iit; 9) collectively as the template located at 9 and qo, q^/2 themselves as the "quadrature" components of the 
template. 

It can be shown, using the stationary phase Fourier transform given in Eq. (|l8|), that 

No = K /2 , (27) 

(«b,? w / 2 > = 0. (28) 

In general, if 

A/# = \\h(t; A = 1, , (29) 

then it can be shown using Eq. ( [I8| ) that 

N$=N = AC/2 ■ (30) 



Note that A = 1 in the above. Eqs. ([28|), (|30| ) also hold to a very good approximation when the numerically computed 
Fourier transforms of the wave forms are used (the typical variation in A/$ is < 1% over $ = [0, 2tt]). By this we 
mean that first the wave forms are generated in the time domain using Eq. (^) and then the Fourier transforms are 
computed using an FFT. Henceforth, we take N<s> to be independent of <& and denote it simply by N. However, N 
does depend on the chirp times via the proportionality constant in Eq. (|l8l). 

At this point it is convenient to define a quantity which we call the strength S of a signal |lz 



S = AN. (31) 

Henceforth, we use S instead of A to parametrise a signal. This quantity is essentially the same as the signal to noise 
ratio [S/N] as defined in ||. 

As a consequence of Eq. (E7h and Eq. (Eq), the test statistic in Eq. (na) reduces to 



A = max 



Cl(x;9') + Cl /2 (x;9')\ 1/2 , (32) 



C (x; 9') = (x(t), q (t-t a ; 9)) , (33) 
C w/2 (a;; 9') = (x(t), q n/2 (t-t a ; 9)) . (34) 

The square root in Eq. (^) is, strictly speaking, not necessary but we retain it in order to make our analysis conform 
to some of the existing literature. 

It has not been possible so far to proceed further analytically with the maximization involved in Eq. (^2|). Although 
several numerical methods are available for the maximization of functions ]l3[| , such methods tend to fail when the 
function involved has many local maxima (as would be the case here), since most of these methods have a tendency to 
converge on one of the local maxima rather than the required global maximum. Also, given the very small expected 
event rate, it is necessary to have a good a priori estimate of the false alarm and detection probabilities and these 
quantities are not easy to compute for such methods. The only alternative left is a search for the maximum over a grid 
of points in parameter space. This method is also simple enough that its performance can be analysed theoretically to 
a large extent. The practical implementation of such a method, called a one-step search, can be motivated as follows. 

For a fixed 9 = 6q : the computation of Co (a;; t a ,9o) (or C n / 2 (x; t a ,6o)), as a function of t a , is equivalent to taking 
the linear correlation JlJ] of x(t) with the template qo(t; 6q) (or q^j 2 in the case of C n / 2 ). Since correlations can be 
computed efficiently using the Fast Fourier Transform (FFT) |]l5| , the maximization over t a alone is straightforward: 
The detector output would be sampled at a rate greater than or equal to the Nyquist rate (in this case ~ 2000 
Hz) yielding a discrete time series x = {xo,x\, . . . ajjv-i ). This time series can then be correlated, using FFTs, with 
analogous time series q and q^/ 2 for the two quadrature components. However, the whole output of such a correlation 
cannot be used since the use of a FFT will yield a circular correlation instead of a linear one. It is, therefore, necessary 
to have a padding of zeroes at the end of each template time series. Let the duration of this padded part be Tp sec for 
some template. Then, for that template, only the first Tp sec of the correlation outputs will be the result of a linear 
correlation and the rest would have to be discarded. Note that Tp will depend on the template parameters since the 
duration of the wave form is parameter dependent (approximately equal to ro + t\ — T1.5). Let the longest duration 
among the template wave forms be £ ma x sec (not be confused with the Newtonian chirp time used in MD96). Then, 
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the shortest linear correlation will have a duration of Tp = T — £ max , where T is the duration of the of the input 
time series x. We will assume in this paper that only the first Tp sec will be kept in each correlation output even 
if the duration of a template wave form is much less than £ max - This appears as a wasteful procedure and a better 
use of computational resources may be possible. However, we do not investigate this issue here since it is not directly 
relevant to this paper. 

Having obtained the correlations with q and q v / 2 i the first Tp sec of each correlation should be squared, corre- 
sponding samples of the two outputs should be added and the square root taken to yield a single time series. We 
call this final time series obtained for some template parameters 9 as the rectified output || of the template 9. Such 
rectified outputs can be constructed for several points in the (ri.s, to) space and the maximum found over each of 
them. Finally, the maximum of all these maxima will yield an approximation to the test statistic A. This, essentially, 
is the scheme of a one step search. We call the set of points in (ri.g, To) space for which rectified outputs are generated 
as the bank of templates. The coordinates of any rectified output sample is given by 9' while that of a template is 
given by 9. 



III. DISTRIBUTIONS OF THE TEST STATISTIC 



To quantify the performance of the detection strategy described above, we need to calculate the probability of a 
false alarm as well as that of detection for a given signal. The former is defined as the probability of the test statistic 
crossing the threshold when a signal is absent in x(t). The latter is the probability of the test statistic crossing the 
threshold when a signal is present. It should be noted that the detection probability need not be the same for all 
signals in a signal family. For instance, large amplitude signals should, obviously, have a larger detection probability 
than weaker ones. We denote the detection probability of a signal with parameters by Qd{Wi ©) where r\ denotes 
the threshold. We denote the false alarm probability by Qo(r]). If the probability density function of A, the test 
statistic, be po(A) in the absence of a signal and pi(A; Q) in the presence of a signal, then, 



Qo(v) 



Hv; ©) 



Po(A)dA , 



pi(A; ®)dA 



(35) 



(36) 



In order to construct po and p±, we start at the lowest level, namely, the density functions of a single sample of a 
rectified output. 

Let z{9') be a sample of some rectified output. Now, this sample will be of the form, z(9') = [C 2 (d') + C^y 2 (#') ] 1//2 , 
where the dependence on x(t) has not been shown explicitly. As mentioned before, both Cq(9') and C^/2(9') are 
obtained from correlations involving the detector output. Thus, they are linear combinations of the samples of x(t) 
and it follows that their marginal probability density function is a Gaussian since the noise is assumed to be a Gaussian 
random process. In the presence of a signal h(t; S, 9' s , 3> s ), their mean values would be, 

C (9') = (h(t; S, q (t-t a ; 6)) , (37) 
C v/2 {6') = (h(t; S, q^ /2 (t-t a ; 9)} , (38) 

while in the absence of a signal, they will have zero means. With our choice of Afo and in Eqs. (25), (p6|), it 

can be shown that the variance of Cq{6') and C 7! /2{9') is unity. It can also be shown, using Eq. (p8|), that Cq(9') 
and Cn/ 2 (9') are statistically independent of each other. Note that, in general, the covariance of Co(0[) and 0^/2(^2) 
need not vanish for 9[ ^ 9' 2 . Given these properties for Co and C^/2, it follows that the marginal probability density 
function of a rectified output sample is a Rician Ri(z) when a signal is present and a Rayleigh R{z) in the absence of 
a signal, 



Ri{z) — z exp 
R(z) = z exp 



-(z 2 + d 2 ) 



I (zd) , 



(39) 
(40) 



where, 



d\9') = C (9') + C v/ 2(9') 



(41) 
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and Io(x) is the modified Bessel function of the first kind (of order zero). For zd ^> 1, the asymptotic form of Ri(z) 
is, 



(42) 



Thus, for z ~ d, a Rician density goes over into a Gaussian. 

Under the stationary phase approximation of Eq. (p"8|), it can be shown that d{6') is independent of the phase <!> 
of the signal. Again, this is a good approximation for the exact numerical case. We assume henceforth that d(9') is 
independent of the signal phase $. Let 



H(e' p ,e' g ) = [{ qo (t-tZ;e p ),q Q (t-ti;e q )) 2 
+ («b(t-*S; 0p),fcr/2(i-tZ; 9 )) 



2 1 1/2 



(43) 



Since d(6*') is almost independent of the signal it follows from Eqs. (P7|),(|3q) and ( |4l| ) that 

d{e') = SH{e', e' s ) . 



(44) 



The quantity H(8' , 6' ) is related to the ambiguity function fiofl . Note that since t a occurs in the phase of a Fourier 
transform, H(9' p , 0' ) depends on t£ and t q a only through At a = t v a — t q a . 

In order to obtain the distribution of the test statistic A we need to know the joint distribution of, in general, all 
the samples. In the presence of a sufficiently strong signal, however, it can be expected that A will occur almost 
always only among those samples of the rectified outputs which have a high value of d. For a typical number of 
samples in a single rectified output of ~ 10 5 , for instance, the location of A is restricted significantly if some samples 
have d > 5.0. Otherwise, A occurs almost randomly anywhere within the output time series. Thus, to obtain the 
distribution of A in the presence of a signal, we need to consider only a restricted subset of all the samples, namely, 
those for which d ^S> 1. The distribution of A would then be that of the maximum over this subset. We describe a 
general scheme for choosing this subset below but for the present, we assume that it is given. As shown above (see 
Eq. (|U)), each sample of this set would have a marginal distribution that is approximately a Gaussian. It is plausible 
that the joint distribution of such samples can also be approximated by a multivariate Gaussian. This possibility can 
be investigated by computing the moments of the joint distribution and comparing them with those of a multivariate 
Gaussian. 

We can express a Rician variable Z = \J X\ + X\ , where Xi is a Gaussian random variable with mean fj,i, as 



I 1/2 



Z = [Oil + 8X x f + (fi 2 + 5X 2 f] ' 

2( fH 5X 1 + fi 2 SX 2 ) 



1 



6 X 2 + 5X1 



-1 1/2 



i4 + i4 



f4 + f4 



(45) 



where 5X± and 5X 2 are zero mean Gaussian random variables with unit variances. In the above expression and in the 
following, we will follow the customary practice of denoting a random variable by a capital letter while denoting its 
value in a particular realization by the corresponding smaller case letter. The probability that 2(fi±SXi + /i 2 SX 2 ) + 
SX 2 + SX 2 be larger than /Ltf + /Lt| can be obtained easily, 



Pr 



2(m8X 1 +n 2 SX 2 ) +5X1 



6X1 



Mi 



> 1 



1 



2ji 



— — / ddexp 



id 2 (v/l + cos 2 [6l] -cos[6l]) 2 



(46) 



where, d 2 = n\ + fi 2 ,. In our analysis, d ~ 7.0 or larger. For d — 7.0, the above probability is 0.008 and it decreases 
rapidly for higher values. Thus, only a small fraction realizations would be such that their binomial expansion in 
terms of (2(/ii<5Xi + fi 2 SX 2 ) + SX 2 + SX 2 )/(fi 2 + /i 2 .) would be non-convergent. It would be a good approximation to 
neglect such realizations and calculate the moments of Z by expanding the RHS of Eq. ( pf5| ) in a binomial expansion 
and taking the ensemble average for each term. 

The same argument goes through for multivariate moments also except for the fact that the fraction of realizations 
for which all the components of the moment can be expanded binomially will decrease as the number of different 
variates increases. For instance, if the third moment {Z\ — ai)(Z 2 — a2)(Z 3 — a^) is required around some point 
(ai, a 2 , a^) then the fraction of cases for which the above expansion will be invalid, taking d ~ 7 for all of them, would 
be ~ 0.008 x 3. In practice there would be a significant overlap of various cases since the Zi would be statistically 
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correlated and this fraction would actually be less. In any case, for low order moments, this method would still furnish 
a good approximation since the fraction of realizations with non-convergent expansions is still small. For higher values 
of d, the fraction of invalid expansions would go down rapidly bringing higher moments also under the purview of this 
method. 

Now consider the restricted subset of rectified output samples mentioned above. Let this set be {Zi, . . . , Z m } 
and the mean values of the Gaussian components, {(Xu, X12), . ■ . , (X m \, X m 2)}, associated with these samples be 
{(^ii)Mi2)) ■ • • , (^1,^2)}- That is, Z { = [Xf l + Xf 2 ] 1/2 and Xu = fj,n and X i2 = \i i2 . Let the strength of the signal 
be S. Following the argument given above, we can express a moment (about mean values) of the joint distribution of 
{Zi, . . . , Zm\ as, 



E[(Z 1 -d 1 ) p ...(Z m -d m ) s ] =E 



1 >> 



(47) 



where E[ ] denotes an ensemble average. Note that the first term is independent of S. Also, since SXn and 5Xi 2 are 
Gaussian random variables, this term is a moment of a multivariate Gaussian distribution. The remaining terms are 
inversely dependent on S as is shown schematically above. In general, the lowest order correction to the Gaussian 
moment will have an S~ 2 dependence for even moments and an S —1 dependence for odd moments. 

Thus, for a sufficiently high S, the moments of the joint distribution function of {Zi, . . . , Z m } are approximately 
the same as the moments of a multivariate Gaussian distribution. This implies that for S> 1, the joint distribution 
of the set {Z\, . . . , Z m } is given by an m-variate Gaussian distribution. It is not easy to see how small the corrections 
to the Gaussian parts of the moments should be for a given error in the detection probability. However, the above 
argument provides a sufficiently strong motivation to proceed with the multivariate Gaussian approximation to the 
joint distribution of {Z\, . . . , Z m }. In order to construct this distribution, we need only compute the covariance matrix 
for {Zi,...,Z m }, 

Suppose we have two sample Zi and Zj having coordinates 6>- = (t*, t{ 5 , Tq) and 8j — (P a ,rf , 5 ,To) respectively. It 
is easy to show, using the stationary phase Fourier transform, that the covariance matrix of { Xu, Xi 2 , Xj\, Xj 2 } is, 
with the columns of the matrix in the same order, 



(10 r s 
1 -s r 
r -s 1 

\s r 1 / 



(48) 



where, 



r= (q (t-f a , 6i), q (t-ti; Bj)) , 
s = (q (t-ti; 9 % ), q v/2 (t-ti; 6 3 )) 



(49) 
(50) 



Both |r| and \s\ are less than unity. The above form of C is also approximately true when the numerically computed 
Fourier transforms of the templates are used. The covariance of Zi and Zj can now be computed using Eq. ( [47| ) with 
the RHS truncated to the first term. We get, 



1 



didj 



[r(fMifiji + A*i2j"i2) + s(nnfi j2 - fiafijij] 



(51) 



The same kind of calculation also yields, 0^ = Uj = 1. The covariance matrix, for the set {Z\, . . . , Z m } above, can 
now be computed using Eq. (|5l]). We can also express the covariance as — \/r 2 + s 2 x, where 



X 



tan 1 




+ tan 1 




— tan 1 


S" 




.M2_ 








. r . 



(52) 



Note that \/r 2 + s 2 is precisely the quantity H{6' i ,6'j). 

We list here the general expressions for the first three multivariate moments obtained using Eq. (47), up to the 
lowest order correction to the Gaussian part. The algebra involved is tedious but a lot of it was automated using 
the symbolic computation package mathematica. Let the three rectified output samples be Z m = [X 2 



X 2 } 1/2 , 



9 



Z n = [Yf+Yg] 1 / 2 and Z D = [Wf + Wf] 1 / 2 and the covariance matrix for (Xi,X 2 ,Yi,Y 2 , Wi,W 2 ) be (columns ordered 
in the same way), 



The moments are, 



where, 



C = 



/ 1 ri si r 2 s 2 \ 

1 — s\ n — ?"2 

ri — si 1 £ u 

si r± 1 — u t 



r 2 -82 t 
\s 2 r 2 u 
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1 / 



(53) 
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(54) 

(55) 

(56) 

(57) 

(58) 

(59) 
(60) 
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[r(fj,afj,ji - iiii^ji) + s(fj,iifj,j 2 + M*2Mji)] 



(61) 



Other moments up to the third order can be constructed from the above expressions by substituting appropriate 
indices. 

Given a bank of templates and the parameters of a signal, the subset {Zi, . . . , Z m } can be chosen as follows. Let 
the coordinate of the signal be 9' s = (t*, r| 5 , Tq). A set of templates is chosen from the template bank which lie in a 
neighborhood of (rf 5 ,Tq), where the size of this neighborhood is adjustable. In the rectified output of each of these 
templates, the sample with the largest value of d is identified. We denote the coordinate of such a sample by 0' a j, 
where the first index stands for the intrinsic parameters (ri.5 and To) of the template and the second stands for the 
location of the sample within the rectified output of this template. From Eq. (|i4|), it follows that the location of this 
sample for a given a can be found by maximising H(9' a k , 0' s ) over the time of arrival k. For each Q' a 3 , we also choose 
2n neighboring samples in the same rectified output, namely, the samples {0' a j +k }, where —n < k < n and k ^= 0. 
Finally, the set of all these samples, namely, the set {9' a ■} and for each (a,j), the set {9' a - +k ; — n < k < n,k ^ 0}, 
gives us the required subset of rectified output samples. Note that H(9[, 9' 2 ) plays a central role in the determination 
of this subset. In our analysis we find that for most of the cases, n = 1 or keeping only the two nearest neighbours to 
9 a j is a good approximation. The choice of the neighborhood of templates is intimately related to the placement of 
the templates themselves and is the subject of the sections which follow. 

The distribution of the test statistic A, in the presence of a signal, is that of the maximum of the set {Zi, . . . , Z m }. 
The joint distribution of this set was shown above to be well approximated by a multivariate Gaussian when the 
strength of the signal is sufficiently high. An analytical form for the distribution of the maximum of a set of Gaussian 
random variables is known only for the bivariate case. There are some approximate methods ]T^ ] for the calculation 
of the distribution but these are impractical for more then four or five variates. However, the distribution for a larger 
number of variates can be estimated from Monte Carlo simulations. A large number of realizations are generated and 
for each realization, the maximum value is recorded and finally an estimate of of the required distribution is obtained. 

Given the covariance matrix Cx of a multivariate Gaussian random variable X = (X%, . . . , Xn), a realization of 
X can be generated as follows. Let {A^; 1 < i < N} be the set of eigenvalues of Cx- Let £ voc be an N x N matrix 
whose ith column is the eigenvector of Cx corresponding to A.;. If W = (Wi, . . . , Wn) is a zero mean multivariate 
Gaussian with a covariance matrix given by diag(l, !,...,!), then 
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Xi 

_ x N 

If X has a non-zero mean vector then this can be added to the RHS of the above expression. It is easy to understand 
the above expression geometrically. In an N dimensional cartesian space, realizations of W will be distributed in a 
spherically symmetric "cloud" . Multiplying the components of W by turns this spehrical distribution into an 
ellipsoidal one. This is the distribution expected for realizations of X in the principal axes frame of Cx- Finally, a 
rotation from the principal axes frame to the actual frame is applied. Another method that can be used |0J is to 
perform the Cholesky decomposition of Cx with the elements of one of the factors chosen in such a way as to give 
the correct covariances for the components of X. We use the method of Eq. (^2|) in our analysis. 

An estimate of the distribution of A can also be obtained using the kind of Monte Carlo simulation that is con- 
ventionally used for studies of parameter estimation accuracy. In such a method, a number of realizations of a noisy 
detector output time series are generated. For each such data segment, rectified outputs are generated for a set of 
templates and the location of the maximum over the outputs is recorded. The distributions of the coordinates of 
the maximum then give an estimate of parameter estimation accuracy. One can also record the values of the max- 
imum and, thus, obtain the distribution of A. Note that a distribution obtained in this way would be free of any 
approximations. 

However, there are some limitations to this method. The first is computational. In a typical simulation in our 
context, each realization of noisy data would have ~ 10 4 samples (for a 1.4, 1.4M binary), for the case of initial 
LIGO, and it would be processed through ~ 5 templates. This leads to ~ 10 6 floating point operations (flop) for 
each realization [fl8|| . A simulation with ~ 2000 realizations would thus involve performing ~2x 10 9 flop (we have 
neglected the cost of generating the noise realization itself). This is not a large requirement computationally but 
when the same calculations are repeated for the advanced LIGO case, for signals having comparable masses, this 
requirement becomes ~ 10 12 flop. This is because the duration of a signal with the above masses is tq ~ 10 3 sec 
for the case of advanced LIGO. Even on a 300 Mflops machine (where an Mflop is 10 6 flop and "flops" stands for 
flop/sec), a typical high end computing power, it would take ~ 1 hr to complete the simulation. In our analysis, 
detection probabilities would be required for various configurations (~ 10 2 ) of template placement and would, thus, 
be quite impractical to compute using this method. There is also a more fundamental limitation. Pseudo-random 
number generators have, in general, a finite period p9[ . For instance, the basic generator provided in the NAg library 
of numerical routines has a recommended maximum output of 4.0 x 10 8 random numbers. For the advanced LIGO 
case, therefore, it is actually not possible to generate more than ~ 200 realizations. This, of course, would lead to 
very poor statistics. 

On the other hand, the method represented by Eq. ( |62| ) is extremely fast and since it does not depend on the 
signal duration, it is equally applicable to both the initial and advanced LIGO case. Let the number of samples in 
the set {0 a .j} be M. Then the total number of samples in the set {Zi, . . . , Z m } would be to = M x (2n + 1). Given 
the covariance matrix for {Z\ 1 . . . ,Z m }, the number of operations required to obtain a single realization would be 
essentially to 2 . Typically, M ~ 5 and n = 2 which leads to 225 flop per realization, a trivial quantity computationally. 
Also, since for each realization only 15 samples need to be generated, the number of trials can be made as large as 
10 7 ! However, we find approximately 10,000 trials to be sufficient for a convergence of the estimated distribution. 
Though the computational requirements for the simulation itself are small, there is a hidden cost in this method, 
namely, the computation of the /i^, r and s. If one were to employ FFTs for their computation, the method would 
again become time consuming. However, these quantities can be computed quite accurately by using the stationary 
phase Fourier transform also. This way, the computation of the covariance matrix also becomes quite fast. Typically, 
the whole simulation including the generation of the covariance matrix takes a few seconds on a 30 Mflops machine. 
This should be compared with the corresponding numbers obtained above for the conventional method. It should be 
noted that this method may be used for simulations of parameter estimation accuracy also. Further investigations in 
this direction are in progress. 

In Fig. ^, we compare the performance of the method of multivariate Gaussian with the exact one. It can be seen 
that the approximate method becomes better as the strength of the signal is increased. In our analysis, a value of 0.95 
for the detection probability will be taken as fiducial and, as can be seen from Fig. [}| the error in the approximate 
method is negligible. 

We now turn to the calculation of false alarm probability. An exact expression for the distribution of A in this case 
is easily obtained when all the rectified output samples being considered form a statistically independent set. In the 
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presence of statistical correlations between the samples, it appears that an analytical treatment is difficult. However, 
it was found in MD96 that Monte Carlo estimates of false alarm probability, as a function of threshold, could be fit 
almost exactly by a formula obtained by assuming that all rectified output samples were statistically independent 
but there were an effectively lesser number of them. That is, if Qo(v) is the probability that the maximum over N r 
rectified output samples crosses 77 (in the absence of a signal) then, 

Q (v) - 1 " exp \-eN r e- r ' 2/2 } , (63) 



-/Vr-eexp 



V 2 



(for 77 > 1), (64) 



where < e < 1. This fit was assumed to hold for higher values of 77 also though they were beyond the range of the 
simulations since, as mentioned earlier, the period of pseudo-random number generators is finite and, consequently, 
it is not possible to generate enough realizations for the estimation of low probabilities. However, we subsequently 
found that this problem was investigated by Rice |2(J (in 1944), though only for the case of a single time series. The 
formulas obtained in that work can also be interpreted in terms of an effective number of samples but the parameter 
e depends on 77 and is not a constant. Specifically, e approaches unity as the threshold is made higher. Therefore, the 
extrapolation of the fit to Monte Carlo estimates that was made in MD96 is not valid. This does not affect the results 
of MD96 significantly, however, because the quantity required in that analysis was 77 for a given false alarm and it 
was shown to be highly insensitive to e. But it should be noted that this implies that Qo(j]) is affected significantly 
by small errors in 77 and, hence, the threshold should be recomputed once the placement of templates is completed. 

Extension of the derivation of Qo(r]) given by Rice (actually, it was the rate of false alarms that was derived) to a 
random field does not appear to be straightforward. Note that the set of rectified outputs from a given template bank 
can be considered to be the samples of an underlying 3-dimcnsional random field, one of the dimensions being the 
time of arrival and the other two being 71.5 and tq. We present here a qualitative argument which can be extended to 
the case of random fields also and show that the same conclusion regarding the behaviour of e should hold. 

The basis of this argument is the assumption that, in the absence of a signal, any rectified output sample is equally 
likely to be the location of A. This can be expected to be true provided the random field is at least wide sense 
stationary. Since the input noise is assumed to be a stationary random process, it follows that the rectified output of 
any one template will also be stationary. However, the random field can be non-stationary because of non-stationarity 
in T1.5 and tq or it is genuinely stationary but sampled non-uniformly in T1.5, tq. For the random field to be wide-sense 



stationary 21 1, it is required that the correlation of any two samples should depend on only their relative displacement 
and not their locations. 

In the present case, the correlation of any two samples, Z\ and Z2, can be obtained exactly in the absence of a 
signal. The derivation given in Appendix ^ leads to the following expression for the correlation, 

ztf* = 2E [H(6[, e' 2 )} - [1 - H(e[, e' 2 f] k [h(9[, e' 2 )\ , (65) 

where E is the complete elliptic integral of the second kind and K is the complete elliptic integral of the first kind. 
Thus, the correlation depends only on H{9' 1 , 9' 2 ). Therefore, if H(9[, 9' 2 ) were location independent, that is if it were 
dependent only on the difference A9 1 between 9[ and 9' 2l then the rectified output field would be at least wide sense 
stationary. 

As discussed below, H(9[, 9' 2 ) is location independent in the case of Newtonian and post 1 -Newtonian wave forms (9' 
is understood to be a different set of parameters for these wave forms) but not in the case of post ^-Newtonian wave 
form. However, for simplicity, consider the case where H(9[, 9' 2 ) is location independent. We apply our argument to 
a single rectified output first. 

Since we will be using an extremely low false alarm probability in our analysis, the threshold required will be quite 
high (typically, 8 to 9a). Suppose that 77 is so high that the probability of two or more simultaneous crossings of 
77, at widely separated locations in the time series, is almost zero. By a wide separation we mean roughly that the 
locations are not closer than the typical correlation length scale. Of course, it is still possible for samples which are 
highly correlated with the one at which A occurs, to simultaneously cross 77. Thus, if Zi is the sample at which A 
occurs (and crosses 77) , then we have assumed above that the probability of Zj crossing 77 in the same rectified output 
is zero, where (|j — i\)At > L and L is the FWHM (say) of the autocorrelation function. We will further assume that 
77 is sufficiently high so that if one of the neighboring samples of crosses 77 along with Z iy it is almost always either 
Zi+\ or Zi_i. 

To compute the false alarm probability, we need to count the number of times A crosses 77 in some N trials for 
N —>■ 00. Under the above assumptions, the number of favorable cases can be counted approximately as follows. First, 
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the number of times a given sample Zi exceeds 77 is counted. For N sufficiently large, it is just the marginal frequency 
n for that sample, 



n = N x / dz z exp 
It) 

N x exp 



z 
~2 



!L 

2 



(66) 



Out of these n cases, some would be such in which either or Zi + \ also cross r\ simultaneously. The number n of 
such cases in which Zi_\ crosses 77 would be n — n x p(zi-i > r\ \ Zi > 77), where p(A | B) is the probability of event 
A given that event B has occurred. The number of cases in which Zi + \ crosses 77 simultaneously with Zi would also 
be n. Now, each sample has n events associated with it that favor A > r\ but out of these, 2n events are common 
with its immediate neighbours and these common events should be counted only once (recall that, by assumption, 
the number of events with more than two simultaneous crossings is negligible). Thus, the total number of events n/ 
that favor a false alarm is, 



7f 

2 



nf~N r xn + N r x (n — 2n) , 

= N r [l -p(zi-i > r] I Zi > 77)] x iVexp 

where, boundary effects have been neglected. Hence, the false alarm probability is, 

Qo(v) = N r [l ~p{z t _i > 77 1 Zi > 77)] exp 
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(67) 



(68) 



A comparison of the above with Eq. (|64|) explains why the latter expression produces a good fit to Monte Carlo 
estimates but now, it can be seen that e = p(zi-i > 77 1 zi > 77) is not a constant as assumed in MD96 but depends on 



exp 



dudv Pzi.Zi^x (u, v) 



(69) 



where, Pz i ,z i _ 1 (u, v) is given in Eq. ( A10| ). Even for such a simple argument, the above expression for e yields values 
that are of the same order as those obtained by fitting the Monte Carlo estimates. For instance, from Eq. ( |69| ) we 
get e = 0.33 for 77 = 6.0 and r 2 + s 2 = 0.9. The typical value for e that was obtained in MD96 was ~ 0.7, for a single 
rectified output. 

It is also clear that the assumptions made above regarding simultaneous crossings of 77 are not strictly necessary. 
The essential point is that simultaneous crossings reduce the number of favorable events and, since the number of 
such events per sample is identical for all samples (under the assumption of stationarity), this can be expressed as 
an effective reduction in the number of samples themselves. In this way, the extension of the above argument to a 
random field is obvious and leads to the same conclusion, namely, that e — ► 1 for 77 ^> 1. In this paper, therefore, we 
use Eq. d63| ) for the false alarm probability but with e = 1. 

In the following, we will need to estimate the threshold that is required to obtain a given false alarm probability. 
For small values of the false alarm probability, we get from Eq. (64), 
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Let the number of templates in the template bank be Nt- Then, 
in a single rectified output. Therefore, 



(70) 



N r = N x Nt, where N is the number of samples 
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(71) 



This shows that the threshold is very insensitive to the number of templates provided the false alarm is kept low. 
For instance, even if there is a relative error of 50% in estimating Nt, the relative error in 77 would just be ~ 0.8%. 
As far as the detection probability of a signal is concerned, such an error is entirely negligible. This point will be of 
importance later in the paper. 
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IV. PLACEMENT OF TEMPLATES FOR A ONE-STEP SEARCH 



In the previous section, we described a method for the calculation of the detection probability of a given signal. 
This method consists of choosing a small set of samples Z from the rectified outputs of templates which are in some 
neighborhood of the signal. The set Z is supposed to be such that the maximum, over all rectified output samples, 
almost always occurs among the members of Z and this requires that each of them should have a high value of 
d (typically, d > 7.0). Thus, the neighborhood of templates should be such that the maximum of each rectified 
output, in the absence of noise, be sufficiently large. 

This motivates the introduction of a quantity called the intrinsic ambiguity function "H(ti 5 ,Tq; Ti 5 ,Tq) which is 
defined as, 

H(6 a , 9 b ) = max H(9' a , 6' b ) . (72) 
t»-t| 

In other words, this is the maximum value that the rectified output of a template 9 a will have if the input consists of 
only a signal 8b having 5 = 1. The role of templates and signals is, of course, interchangeable here. For 5^1, the 
maximum value will simply be STi.. We term this value as the observed strength S Q b s of the signal. 

Clearly, the larger the "width" of the intrinsic ambiguity function, the more sparsely can templates be placed around 
a signal in order to obtain the same detection probability. In this sense, the intrinsic ambiguity plays a central role in 
the determination of the density of templates and thus the computational cost of a one-step search. We first discuss, 
in the following section, the calculation of TL and some of its relevant properties. This is followed by a discussion of 
template placement for a one-step search. 

A. The intrinsic ambiguity function 

Using the stationary phase Fourier transform given in Eq. fllq) , it can be shown easily that, 

H(9' a , 9' b ) = i [FLKJ' b ) + FL(0a,8' b )] 1/2 , (73) 

F cos (e' a ,9' b ) = ^ s^r 7/3 C0S M/)] . ( 74 ) 

F sin (9' a ,9' b ) = J" -^y/~ 7/3 sin[«(/)] , (75) 
«(/) = 2ir[At a f + A7*Vk(/)] ; k G [0, 1, 1.5], (76) 

k 

' 3 =r*f7> r7,s ' <77) 

where At a = t a a — t h a , At]~ = r£ — r|. The quantities F cos /f3 and F s i n //3 are nothing but the quantities r and s defined 
in Eq. (|49| ) and Eq. ( |50|) but expressed in terms of the stationary phase Fourier transform. They can, of course, be 
calculated exactly by generating the wave forms in time domain using Eq. @ and computing their correlations using 
FFT. We find that F cos and F s - ln reproduce the corresponding exact quantities quite faithfully for both the initial 
and advanced LIGO noise spectral densities. This also holds to a large extent when f c is replaced by the least of 
the plunge cutoff frequencies corresponding to the two wave forms. Thus, H(9' a , 8' b ) calculated using Eq. ([73]) also 
reproduces faithfully the corresponding exact results and this would also be true for H. provided the location of the 
maximum in Eq. (|72|) is obtained accurately. 

We obtain the location of the maximum in Eq. ( |72| ) in two steps. First, an initial estimate for the required value of 
At a is obtained as described below. This is followed by a search for the true maximum around this initial guess using 
a bracketing and golden search algorithm (MNBRAK and GOLDEN in |k|). In order to get the initial estimate, we 
solve the integrals in Eq. (|74|) and Eq. ( f75| ) using a stationary phase approximation but with the point of stationarity 
chosen in such a way as to yield the maximum value for the RHS in Eq. (|7^). Let the desired stationary point be 
/ = fo- Then, for fixed values of Ati.s and A to, the condition of stationarity yields 

A* o = -(At + Ati-Ati.b) 

+At x 8 + At iX 6 - An. 5 x 5 , (78) 
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where, x — [/o//a]~ 1/ ' 3 - Substituting the above back into the integrands in Eq. (|7^) and maximising the resulting 
expression over fo yields the required value. This value of /o is obtained only once for a particular choice of Ati.5 
and Ar . For any other (Ati. 5 , Ato), the same value of / is used to obtain an initial estimate for At a using Eq. (|7g|). 
The algorithm given above is quite fast as compared to the exact calculation. 

Contour plots of H(9 a , 9b), as a function of Ob with 9 a fixed, are shown in Fig. [|. Also shown are the eigenvectors 



of the Hessian of H which is defined as, 

H 



1 d 2 H(9 a , 9 b ) 



i 'j » 



2 391891 



(79) 



Since H(6 a , 9 a + A9) is maximum at A9 = 0, the 7i surface is quadratic in a sufficiently small neighborhood of 9 a 
and, as shown in Fig. ^, the inner-most contours are elliptical. The orientation and axes lengths of such an elliptical 
contour can be obtained in terms of the eigenvalues and eigenvectors of the Hessian. Let \i(0 a ), \2{9 a ) be the smaller 
and larger eigenvalues of H(0 O ) respectively and let e\ a and e^a be the corresponding eigenvectors (normalized to 
unity). Then, the length of the semi- minor and semi-major axes of a contour at level e are given by, 



h(e,9 a ) = ^^=, (80) 

\/-M{t>a) 

WMJ = 7& (81> 

respectively while their orientation is given by the respective eigenvectors. Note that since the eigenvalues are negative, 
|Ai(0 a )| > |A2(0 a )|- Since, we already have a nearly accurate method for computing H, Hij can be calculated simply 
by approximating the derivatives in Eq. (f7^) by finite differences. It should be noted that at a given location, the 
directions in which the contours at successively lower levels are most elongated suffer progressively larger rotations 
with respect to these eigenvectors. We call this the shear of the contours. Typically, Xi(9 a ) and ej a provide a good 
estimate of the size and orientation of the contours for Tt > 0.95. 

The intrinsic ambiguity Ti. is not independent of its location in parameter space. That is, if A9 = (Ati.5, Ato) be a 
displacement vector, Tt(0 a , 9 a + AO) ^ T~t(0b, b + AO) in general for a ^ b . For our purpose, the most appropriate 
way of characterizing this location dependence would be to investigate the change in the dimensions of the inner-most 
contours of H. This is because we are primarily interested in the detection of signals with low strengths and for such 
signals the templates in a one step search would be placed closely. For instance, we find from earlier works ( MD96 
and ||]) that for low values of S, the spacing of templates is such that the signal with the lowest detection probability 
has U(0 t , 9 S ) ~ 0.97, for some template 9 t . In Figs. §, § and Figs. |, |, we have plotted Zi(0.97, 9 a ), Z 2 (0.97, 9 a ), the 
area of the ellipse irh I2 and the angle between the T1.5 axis and e\ a as functions of a - The lowest contour level in 
each plot is close to the minimum value which that quantity takes over the whole of the space of interest. 



B. The geometrical configuration of a one-step template bank 

In MD96, we had introduced a set of criteria which a template bank for a one-step search was required to satisfy. 
These criteria were (CI) every signal, having a strength S greater than a given minimum strength S'min, should have 
a detection probability greater than a given minimum detection probability Qd,min and (C2) The false alarm should 
stay below a specified level Qo,max- Throughout the following, Qd,min = 0.95 and Qo,max is chosen to be such that 
the average rate of false events is 1 event/y. Apart from the above criteria, we also demand (C3) that the templates 
be spaced as sparsely as possible so as to minimise the computational cost. 

An obvious way to fix the template placement would be to search through all the possible placement configurations 
and find the one satisfying the three criteria stated above. We have already introduced formulas for detection and 



false alarm probabilities in Section III which can be used in checking CI and C2. However, such a blind search in 
configuration space is impossible to perform in practice since the number of templates can be expected to be quite 
large. Instead one can make some reasonable assumptions regarding the geometry of the final configuration and 
then proceed to perform a limited search within that framework. We will now present an argument that suggests an 
approximation to the optimum geometry of template placement in the post ^-Newtonian case. This approximation 
should be good enough for estimating the performance of a two-step hierarchical search but a more careful analysis 
would be required when such a scheme is actually implemented. 
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1. The case of location independent TL 



Consider, first, a simple hypothetical situation in which the intrinsic ambiguity is location independent, that is, 
H(9 a , 9 a + AO) = H(A6). Recall that the detection probability of a signal was determined by a subset Z of samples 
belonging to the rectified outputs of templates in some neighborhood of the signal. Thus, the detection probability of 
a signal is determined entirely by the local distribution of templates around it. Let the coordinates of the signal be 9 S 
and that of the templates be {9 S + A6%, . . . , 9 S + A9 p }, where P is the number of templates in the neighborhood. Also, 
the samples in Z which contribute most to the detection probability are the ones which are located at the maxima of 
the rectified outputs (for zero noise) of these templates. Let the set of such samples be Z' = {Z^} C Z, i — 1, . . . , P. 
The remaining samples typically contribute only a few percent more to the detection probability. Then, 

Z T l ~STL(9 s ,9 s + A9 l ). (82) 

Now, if TL is assumed to be location independent, two different signals (having the same strength) would have the 
same detection probability provided the local configuration of templates {A6>i, . . . , A9 P } around them is the same. 
Strictly speaking, this statement is not true unless H(9' a , 9' b ) is also independent of location. This is because the 
covariance matrix (determined by H) of the set of samples Z' need not be location independent even though their 
mean values (given by TL) may be so. However, we will proceed with the assumption that variation in the covariance 
matrix is a negligible effect. The validity of this assumption can be checked after the final results have been obtained, 
as we do below. 

An intrinsic ambiguity function which is location independent is not unrealistic since the Newtonian ||[l2| as well as 
the post^Newtonian Q wave forms are known to have such an intrinsic ambiguity for the right choice of parameters. 
In fact, in these cases, the function H(8[, 9' 2 ) itself is dependent only on A9. Thus, the detection probability of a signal 
in the case of Newtonian or post ^Newtonian wave forms depends strictly on the local configuration {Af?i, . . . , A9 P } 
alone. 

Now consider a configuration of templates in which the distribution of templates in the space of interest is inhomo- 
geneous but which is claimed to satisfy CI, C2 and C3 for some >!? m m- It then follows that the detection probability 
of a signal with S — S mm in a region sparsely populated by templates is at least Qd,mm- However, this implies 
that a region where templates are spaced densely is over populated because, as shown above, the same sparse local 
distribution should suffice everywhere. Therefore, a further reduction in computational cost can be brought about 
by removing some of the templates from the over populated region. This implies that an inhomogeneous distribution 
cannot be an optimal one. Hence, when TL is location independent, templates should be distributed homogeneously 
which is equivalent to placing them on a regular grid. 

A two dimensional regular grid is specified by a single unit cell. In the present paper we will assume that this unit 
cell is a parallelogram. Such a unit cell is specified by the lengths of two adjacent sides, h and I2, and the angles, 
ail, Q!2, that l\ and I2 make with some reference direction (see Fig. ^|). To find the optimal placement of templates, 
therefore, a search can be performed in (Zi, I2, ot\, 012) space for the unit cell having the largest area under the constraint 
that the resulting grid of templates satisfy CI and C2. Such a search would be computationally expensive even with 
the fast methods that were introduced in the earlier sections. Note, however, that the unit cell with the largest area, 
for given li and I2, is a rectangle. Hence, for given li and I2, a search should first be performed among all rectangular 
unit cells to locate the ones that satisfy CI and C2. This is equivalent to just a rotation of the grid which involves 
only one of the angles and, hence, saves significantly on computations. Once the largest rectangular unit cell that 
satisfies CI and C2 has been found, the search can then be extended to non-rectangular unit cells with larger areas. 

The computational cost can be reduced further by making an educated guess for the orientation of the rectangular 
unit cell. For instance, let the contours of TL be ellipsoids with the same orientation. Then it can be seen heuristically 
that the largest rectangular unit cell should be obtained when l\ and I2 are oriented along the major and minor axes, 
that is, the eigenvectors of the Hessian H. If the contours exhibit a shear, then the largest rectangular unit cells for 
different values of S^un would, of course, be oriented differently. In such a case also, the computation involved in the 
search can be reduced by starting with an orientation given by the eigenvectors of H and then searching a small range 
of ot\ around this orientation. 

In order to check that the above argument is reasonable, consider Fig. [l0| For a rectangular unit cell oriented along 
the T1.5 and To axes and having arbitrary dimensions l\ and I2, we show the detection probability of signals, having the 
same strength, which lie in the interior of the cell. The threshold has also been chosen arbitrarily and samples for the 
set Z were chosen from the rectified outputs of the four templates at the vertices of the unit cell. Also superimposed 
on this detection probability map are the contours of TL for one of the templates. It is clear from the figure that, as 
expected, the detection probability map closely follows the contours of the intrinsic ambiguity. Roughly speaking, 
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the detection probability contours are formed by the "overlap" of the Ti contours. Therefore, it can be expected that 
if the unit cell were oriented along the eigenvectors of H, then the area can be made larger, keeping the minimum 
detection probability the same, because this orientation would maximise the overlap of the Ti contours. In Fig. [h] 
we also show the detection probability map, for the same threshold as above, when the unit cell is oriented along 
the eigenvectors of H (l\ along e\ a and l 2 along e 2a ) . The minimum detection probability is much larger now which 
implies that the area can now be increased further. In the following, we will restrict the parameter space for the unit 
cell to be only (1%, l 2 ) and orient the sides along the eigenvectors of the Hessian. That is, the sides of the unit cell are 
given by h e la and l 2 e 2a . 



2. The case of post 1 ' 5 -Newtonian Ti 

Consider the family of post 1 ^-Newtonian wave forms now. As shown in the previous section, Ti is not location 
independent in this case. However, if this variation is small over the scale of few unit cells then it can be expected 
that the optimum template placement for the post ^-Newtonian case would also be close to a regular grid. In earlier 
works (MD96 and Q), the typical spacing of templates for low values of S m i n turned out to be such that for any 
signal, the value of Ti was at least ~ 0.97 in some template. In the present case, if templates are placed along the 
eigenvectors of the Hessian, this would imply that the typical lengths for the sides of the unit cell be ~ 2 Zi(0.97, 9 a ) 
and ~ 2?2(0.97, 9 a ). Thus, the effect of the location dependence of Ti on the placement of templates can be studied 
by comparing the change in li, over a distance 2^ along e*i a (that is, (0.97, 9 a + e"j a ) — (0.97, a )), with the value 
of li at that point. For the scales used in Fig. || and Fig. |t], a line along any of the two eigenvectors would be nearly 
horizontal at any point. Therefore, the change in li along ej a can be obtained approximately by simply measuring 
the change over a line of constant tq. In Figs. O, |l2|, we plot the quantity 



M0.97,(t°,ti, 5 + 2Z?)] 
y0.97, (r °, n.5)] 



100 , (83) 



where I® is the value of k at the intersection of the tq = Tq line and the boundary AB or BC as the case may be [E2|. 
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A similar quantity 5 a can be constructed for the angle a\ that the semi-minor axis makes with the T1.5 axis 



5a = 



x 100 . (84) 



ai(To, n.5) 

This has also been plotted in Figs, [ll], [l^. 

It is evident that the variation of the dimensions and orientations of the unit cells, over the scale of a single unit 
cell itself, is quite small over a large portion of the space of interest for both initial and advanced LIGO. In general, 
this variation becomes more rapid towards the high mass, or low to, region. However, it is still small for advanced 
LIGO though it may have some significant effect in the case of initial LIGO. Thus, at least in the case of advanced 
LIGO, one can expect that the optimum placement of templates will be along an "adiabatically" changing grid over 
most of the space of interest. Since the rate of detectable events is not expected to be large for initial LIGO, we 
will not investigate the placement of templates for initial LIGO any further here. Instead we will concentrate on 
advanced LIGO and assume that a quasi-regular grid will be obtained for initial LIGO also. Further, in the following 
analysis, we will approximate the quasi-regular grid above by a set of piecewise regular grids. That is, a set of patches 
covering the whole of the space of interest where the unit cells in each patch are identical but differ in dimension and 
orientation in different patches. 

Though, in principle, the unit cell in each patch can be determined by using the algorithm given earlier, this would 
again be impractical because now a placement algorithm would have to search a two dimensional parameter space, l\ 
and I2, for each patch. However, if the assumption that the detection probability is almost completely determined by 
Ti alone were true, then the search would collapse to just two dimensions. This can be seen as follows (we call this 
assumption Assl for convenience in the following). For a small displacement A9 = x± e± a + x 2 e.% a at 9 a , 

Ti(9 a , 9 a + A9)~l- [Xi(9 a )xl + X 2 (9 a )x 2 2 ] , (85) 

where we have rotated the coordinate system locally so that H(0 O ) is diagonalized. Consider a different location 
9b where Aj(#b) = ai\i(9 a ). Then, for a small displacement at 9^ (in a rotated coordinate system that diagonalizes 
H((9f,)) and at the same order of approximation as in Eq. (pq) , 

H[9 b , 9b + J2 x £ib/\fa~i] = "H[9 a , 0a + J2 x &*\ ■ ( 86 ) 
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Now, suppose that the unit cell at 9 a , having dimensions l\ and I2, satisfies CI for a given threshold. This implies that 
every signal in the interior of the cell has a detection probability greater than Qd.min- Let the relative coordinates of 
such a signal be 9 S — (ei h, 62 h), where < e» < 1. If Assl is true, then the detection probability must depend only on 
Hi = H(6 a +6 S , 9 a ), H 2 = H{9 a +9 S , 6 a +he la ), H 3 = H(6 a +6 S , 9 a +l 2 e la ) andft 4 - H{9 a +9 S , e a +he la +l 2 e 2a ) (we 
have neglected the contribution from other templates for the present but this does not affect the argument). It then 
follows from eq. (|8^) that, for the same threshold, the detection probability of a signal with relative coordinates 
(ei l\l \fai, €2 I21 1 \[c*2) at 6*b will be the same as that of 9 S at 9 a . Hence, for a given threshold, if a unit cell with sides 
h, I2 at 9 a satisfies CI, then so would a unit cell with sides h/y/cn, h/y^l at 9b- 

Assume that for a given threshold, the largest unit cell that is compatible with CI is unique (note that the orientation 
of unit cells has already been fixed and only rectangular unit cells are being considered). We call this assumption 
Ass2. Now, suppose that the optimum solution compatible with all the three criteria CI, C2 and C3 has been 
obtained by some means. That is, the sides of the largest unit cell compatible with CI in each patch as well as a 
common threshold have been found. From Ass2 and Assl, it then follows that if the unit cell dimensions are l\ 
and I2 in any one patch, the dimensions of a unit cell in any other patch must be l\j \fcl\ and I2I \fot2- Hence, when 
searching the 2N p parameter space of unit cell dimensions (N p being the number of patches), only the subspace l\, 
I2 for any one unit cell needs to be searched. We should emphasize here that the above argument is by no means 
a rigorous proof. Given the complicated interdependencies of various quantities (for instance, even the number of 
patches and also the extent of a patch may depend on the dimensions of the unit cells), it would be difficult to cast 
the problem into a tractable mathematical form. However, we find the above argument sufficiently suggestive and 
the conclusions reached as plausible. The assumption Ass2 is actually not required since it can be expected that the 
total number of templates, hence the threshold for a given false alarm, will only depend on the areas of the unit cells 
in each patch and not on their individual dimensions. Thus, one can always choose the optimum solution to be the 
one where unit cells are scaled verisons of each other, without violating C2 or C3. 

How large can l\ and l 2 be in the post 1 ' 5 -Newtonian case before the detection probability for scaled unit cells starts 
showing significant errors? We have checked this empirically and the results are presented in Figs. |l3| (initial LIGO) 
and Figs. |lj (advanced LIGO). In each figure we present our results for different values of 5 m in which are chosen 
to encompass the typical range of S m i n that will be considered later. The detection probability that we consider 
throughout our analysis is 0.95. Hence, we compute the errors in detection probability at the 0.95 level. Also, it is 
not enough to compute the error for only one signal point since it may depend on the signal location. Therefore, the 
maximum error among three different signals is shown. 

For each plot, we take two widely separated locations 9a and &b- The x-axis and y-axis are the values of l\ and 
I2 at 9a- The corresponding quantities at 9b are l[ = h/y/jx and 1' 2 — I2/-J12, where 7^ = \i{9 b) / ^i{9 a) ■ At each 
location the detection probabilities of three representative signals are obtained (this anticipates the discussion of the 



Section gVTj), namely, the signals 9i a = 9 a + {h ei a + l 2 e 2a )/2, 9 2a = 9 a + liei a /2 and 9 3a = 9 a + l 2 e 2a /2, where a = A 
or B. Let the threshold at which the detection probability of 9iA equals 0.95 be r/i. Let the detection probability of 
9iB at the same threshold rji be Qdi- The quantity plotted on the z-axis is the maximum relative error among the 
three signals. That is, 

Qdi 



max 



1 - 



0.95 



x 100 



Thus, we are plotting the maximum relative error in the detection probability (at the 0.95 level) as a function of the 
unit cell dimensions. As mentioned earlier, the typical one-step spacings that can be expected are 2Zi(0.97, 9 a ) and 
2 ^(0. 97, 9 a ). For these values, we see from the figures that the typical error is < 2%. In fact, the errors stay small 
for much larger values of the unit cell dimensions. Hence, for a one-step template placement involving low values of 
S'min, Assl can be assumed to be valid for post 15 -Newtonian wave forms. 



C. The number of templates for a one-step search for post 1 5 -Newtonian wave forms 

In order to obtain the computational cost of a one-step search as well as the threshold for a given false alarm, 
the number of templates in a grid has to be obtained. This is not a straightforward task, however, because of the 
non-trivial shape of the boundary and the variation in the area of the contours with change in location. 

The non-trivial shape of the boundary would cause some templates in any regular grid to fall outside the space of 
interest. However, as can be seen from Figs. (|l|), (||), in which the unit cells are almost horizontal because of the scales 
used, such an effect would be significant only near the region around vertex A. In this region, even a single unit cell 



18 



may span both the boundary segments AC and AB. Note, however, that the segments AB, as well as BC, are not 
strict limits. That is, astrophysically valid templates can also exist beyond them. This is not true, however, for the 
segment AC which is the image of the principal diagonal in the (mi, 7B2) plane. Thus, although no template should 
be placed on the left of AC, it is acceptable if some templates in a grid breach AB or BC. This allows a quasi-regular 
grid to be placed in the region near A, as shown schematically in Fig. [l~5|. This patch of templates can be used to 
cover the space of interest from A till that value of tq at which the width of the space of interest becomes comparable 
to the length I2 of the unit cell. The number of extra templates thus added will not be significant compared to the 
total number of templates that will be required to cover all of the space of interest. 

The effect of the variation in the area of the TL contours on the number of templates can be incorporated approx- 
imately as follows. Recall that in the previous section we showed that unit cells in different patches can be taken 
as scaled versions of a standard unit cell, where the scale factors were the ratio of corresponding eigenvalues of the 
Hessian. Thus, the area of a unit cell will vary with location according to these scale factors only and, hence, the 
relative change in the area will be the same as that in the area of the 0.97 contour which is shown in Fig. ^ and Fig. @. 

Let l\ and I2 be the dimensions of the unit cell in that region of the space of interest where the variation in the area 
is small (say, maximum relative change of ~ 15 %). For initial LIGO this is roughly the region between vertex A and 
the contour at 0.002 sec 2 in Fig. [| while it is the area between A and the contour at 0.04 sec 2 for the case of advanced 
LIGO. Let Ac 1 -c 2 be the area between contours c\ and C2 and a = l\ li- Then, for the case of initial LIGO, the 
number of templates that lie in Aq. 002-0. 003 would be approximately Ao. 002-0. 003/(1-5 a) since the area of the contour 
increases by ~ 50 % in this region. Similarly, the number of templates in Aq, 003-0. 004 would be ~ Ao. 003-0. 004/(2-0 a) 
and so on. Let Ac?i-C 2 — Ad-Cs/A, where A is the area of the whole of the space of interest (see Eq. jl4|)) and Nj, 
be the number of templates in the region where the variation in area is fast (i.e., below the 0.002 sec 2 contour). Then, 
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where we have not taken more terms because their corresponding areas are negligible (even Ao. 007-0. 008 = 
0.07^0.002-0.007)- The values of A?i-c 2 are, Ai. 002-0. 003 = 0.193, A. 003-0. 004 = 0.090, A. 004-0. 005 = 0.061, 
A. 005-0. 006 = 0.048, A). 006-0. 007 = 0.041 and Ao. 007-0. 008 = 0.031. We call the coefficient of A/a on the RHS of 
Eq. ( |87j ) as k. For the case of initial LIGO, therefore, k = 0.233. Similarly, for the case of advanced LIGO, 

A). 04-0.05 . A. 05-0. 06 . A. 06-0. 07 
K = , 

1.25 1.5 1.75 ' 

where Ad. 04-0.05 = 0.183, A).05-o.06 = 0.062 and A.06-0.07 = 0.020 which give, K — 0.199. To understand what these 
values for k mean, assume that the number of templates in the remaining region of the space of interest (that is, the 
region with a slow variation) can be obtained by simply dividing its area by that of the unit cell. Then the total 
number of templates iV^ would be, 

i\£ = (89) 

For initial LIGO = °- 464 > which S ives n t - 0.769 A/a. Similarly for advanced LIGO, N f T ~ 0.934 A/a. This 

clearly shows that the variation in the area of unit cells has a small effect in the case of advanced LIGO. 

We combine the two approximations discussed above to give the following algorithm for estimating the total number 
of templates. Recall that in the region of large To, all the three quantities l\, I2 and a.\ vary quite slowly. For the 
purpose of counting the number of templates, therefore, we will assume the orientation and dimensions of a unit cell 
to be constants in the regions (a) between vertex A and the 0.002 sec 2 contour for initial LIGO and (b) between 
vertex A and the 0.04 sec 2 contour for advanced LIGO. We choose an average value of ot\ = 38° for initial LIGO and 
ct\ = 45° for advanced LIGO, where a.\ is the angle between the semi-minor axis and the T1.5 axis. The final results 
are quite insensitive to the choice of these angular values. In the first step of the algorithm, we count the number of 
templates in the region near A by placing unit cells as shown schematically in Fig. [l5| The unit cells are "stacked" 
below each other until the length of the segment along a tq = const, line equals h- Let this value of tq be Tq 9 ^) and 
the number of templates thus obtained be N^ q . The area, A eqi between the vertex A and the To = Tq 9 (/2) line is then 
found. The total number of templates is then obtained as, 



- + K-. (90) 
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The output of this algorithm is shown in Fig. |l6| and Fig.|l7|, where we have also shown the values obtained if the 
number of templates is estimated simply as A/a. Again, it can be seen that the effect of variation in Ti. contours is 
small for advanced LIGO. 

Almost all the problems associated with the non-trivial shape of the boundary of the space of interest can be 
eliminated if instead of the segment AB, a rectangular corner ADB were used. That is, D has the abscissa of B 
and the ordinate of A. However, the region between the segment AB, AD and DB is then mapped, in the (mi, 7712) 
plane, onto a negligibly small area. Thus, although the number of templates will increase significantly, not much will 
be gained in terms of the range of detectable binary systems. 

The boundary can also be made simple by going over to a different set of parameters, like the masses (mi, 7712). But 
we found that in such cases the intrinsic ambiguity function shows excessive location dependence. However, there may 
exist a coordinate system in which both the boundary of the space of interest is simple and the intrinsic ambiguity 
function does not show much variation. This approach needs to be explored more thoroughly. The problems with the 
counting of templates, discussed above, are also present for the coordinates used in @,@|. 

To summarize, the following conclusions emerge from the discussion presented so far. (i) In a case where the 
intrinsic ambiguity is location independent (as happens, for instance, in the Newtonian and post^Newtonian case), 
templates should be placed on a regular grid (neglecting boundary effects) . (ii) The unit cell of the grid should have 
the largest area (in order to satisfy C3) while satisfying CI and C2. Assuming that the unit cell is a parallelogram, 
we gave a practical algorithm to find the four parameters (li, 1%, a\, a.2) of this optimum unit cell, (iii) The location 
dependence of Ti. in the post 15 -Newtonian case is quite weak over most of the space of interest (at least for the case 
of advanced LIGO) as shown in Fig. [ll]and Fig. |l~2|. This implies that the placement of templates in this case should 
also be on an approximately regular grid, (iv) We can make a piecewise approximation to this grid where each piece, 
or patch, is regular (formed by translating the same unit cell), (v) If the detection probability of a signal were to be 
determined almost completely by the intrinsic ambiguity, then only a single unit cell for any one patch needs to be 
determined. We checked that this assumption is true for the post ^-Newtonian wave form, (vi) Since the segment 
AB is not a strict boundary, it is acceptable if some templates from a grid near vertex A breach it. We, therefore, 
put a quasi-regular grid of single unit cells "stacked" vertically near A. In the remaining region, variation in the area 
of unit cells was approximately taken into account while estimating the number of unit cells. We emphasize here that 
these conclusions will not hold for large values of S m i n but only for values that are sufficiently low so as to make the 
unit cells small. 



D. Algorithm for the determination of the optimum unit cell 

We now present the algorithm that we have used in this paper for the determination of the parameters l\ and I2 of 
the optimum one-step unit cell. Although the algorithm that was obtained earlier is practical enough, we can simplify 
it further as follows. 

First, as has already been shown earlier (Section |n|), the threshold is very insensitive to the number of templates 
when the required false alarm is small. Since, as shown in Fig. [Hi] and Fig.^, the relative error between the approximate 
count A/ a and the "exact" count is ~ 30% or less, the threshold is affected negligibly if the approximate count is 
used in its determination. 

However, it must be emphasized here that the false alarm is very sensitive to changes in the threshold and, therefore, 
ultimately the threshold should be fixed very accurately. For instance, let the number of templates be 3 x 10 5 and 
the length of a data segment be 8192 sec with a sampling rate of 2048 Hz (typical values for advanced LIGO). Then, 
for an average false alarm rate of 1 event/y, the required threshold is 8.661. If an error of, say, —5% is made in 
the determination of this threshold, then the false alarm rate becomes ~ 38 events/y while if the error made is 
+ 5%, then the false event rate falls to 0.02 events/y. The latter situation is definitely more preferable, however, 
since the detection probability of signals will not drop too much (typically by 5%). This in turn implies that an 
overestimation of the number of templates is better than an underestimation and that is precisely what happens when 
the approximate count of template is used (provided the effect of the boundary near vertex A is negligible, which 
should be so for small unit cells). 

Secondly, given a unit cell, it is sufficient to check that signals that have the least detection probability be detectable 
with a probability > Qd.min = 0.95. For a rectangular cell, there exist three such signals, namely, the signals at the 
mid-points of l\ and I2 and the signal at the the centroid of the rectangle. This is also borne out by Fig. |l^ (lowermost) 
which shows that these three signals lie at the minima of the detection probability map. This is also the reason that 
such a set of signals was used in Figs. OL O. 
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The Algorithm for one-step template placement: 
(i) The value of S m i n is fixed. We make a few preliminary coarse runs to find out the range of values of S n 
the unit cells are sufficiently small (so that the shear is unimportant). As mentioned earlier, we keep Qd.min = 0.95 
and Qo,max is kept such that the average rate of false events is 1 /y. Let the duration of each input data segment be 
T sec. Then, 

Q ,max = {T ~ 6nax)/(365 X 24 X 3600) , (91) 



where £ max is the duration of the longest template (see Section [II B| ). Recall that we denote T — £ max by T P . 

(ii) We choose a point in the (ri.5, To) space such that when a unit cell is constructed around it, all the templates lie 
well within the boundary. For instance, in the case of initial LIGO, we choose the point (1.3,50.0). Let this point be 

(iii) The unit eigenvectors e\ a and eia of H are found at Q a . We consider rectangular unit cells such that (a) 9 a is 
always the same vertex for all of them (b) the sides are l{e\ a and lie^a with l\ > 0, I2 > 0. The values of l\ and 1% 
are chosen to lie on a regular grid in (l\, l 2 ) space. Typically, we keep (a) l\ £ (0.01, 0.05) sec and I2 € (0.05, 0.2) sec 
for the initial LIGO (b) l\ £ (0.04, 0.1) sec and I2 £ (0.2, 0.6) sec for the advanced LIGO. The number of grid points 
is kept at ~ 10 x 10. 

(iv) Given a point (h, 1%), the threshold 77 required to get a false alarm Qo.max is computed using Eq. (|63]). The 
total number of rectified output samples is N r — NX x Tp x 2048, where Tp is the padding in the time series of the 



template with the longest duration (see Section II B) and is the total number of templates which is taken as 



t Area of the space of interest A 

1 V 'T 1 ■ ■ ■ p I J 

area of the unit cell = 

(v) For each unit cell we consider the signal points S\ = 8 a + (h/2) ei a , S2 = a + (h/2) &2 a and S3 = 9 a + (h/2) e% a + 
(?2/2)?2a- Their respective detection probabilities Qd,i, Qd,2 and Qd,3, for the threshold 77, are computed using 



the algorithm based on a multivariate Gaussian joint density (see Section III). The set Z is chosen from among the 
rectified outputs of the templates at 6 a , 9 a +hei a , Oa+h^a, Oa+heia+h^a, 8a+heia—he2a and 8 a ~h ei a + l 2 e 2a . 
The last two templates are included in order to take care of any contribution to the detection probability due to the 
shear in Ti. contours. 

(vi) If the minimum among {Qd.i, Qd,2, Qd,3} is larger than Qd.min then the point (Zi, I2) is recorded or else not. Let 
the set of unit cells which qualify thus be L. 

(vii) The unit cell with the largest area among L is chosen as the optimum unit cell. 

We do not proceed further to find a larger non-rectangular unit cell because for the values of >!5 m in used, there will not 
be much of an improvement. 



V. TWO STEP HIERARCHICAL SEARCH 



Our scheme for a two-step hierarchical search involves the use of two banks of templates, in both of which we use 
the same family of templates. In one of the banks, templates are spaced sparsely in (intrinsic) parameter space and 
this is called the first stage bank. The second stage template bank consists of templates placed more finely. The 
detector output is first processed through the first stage templates and the locations of those templates are noted 
in whose rectified outputs there was at least one crossing of a threshold r\\ (the first stage threshold). The value of 
r/i is kept sufficiently low so that, even though the templates are sparsely spaced, all signals (with strength greater 
than some S min ) can produce at least one crossing in a nearby template with a probability of ~ 0.95. In the next 
step, for each of the first stage templates that produce a crossing of 771, the detector output is processed through a 
neighborhood of second stage templates around it. The maximum over these second stage outputs is then compared 
with a second stage threshold 772 to check for a detection. In this way a significant saving in computation occurs since 
the number of templates used on the whole is much less than if the second stage bank alone were used. 

The first stage templates cannot be spaced too coarsely, however, because 771 has to be lowered and at some point 
the number of crossings because of noise alone becomes large. Since each such false crossing would involve the use of 
second stage templates, the computational cost starts rising if the first stage templates are spaced too coarsely. Thus, 
there is a non-trivial optimization problem that needs to be solved while setting up a two-step hierarchy. 

It was shown in MD96 that the correlations between templates allows a false event (crossing of 772 due to noise 
alone) to slip through the hierarchy in spite of the presence of two thresholds. This is essentially because of the fact 
that the hierarchy is designed to allow the easy passage of a signal and if a noise realization is such as to produce a 
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crossing of 772, which is quite high (typically ~ 8.0), then it would have sufficient "resemblance" to some signal (in 
its phase information) to allow it to pass through the hierarchy. Stated in another way, this is because, for a first 
stage template and its neighborhood of second stage templates, the crossing of ry 2 is not statistically independent of a 
crossing of 771 . This implies that the second stage template bank and threshold should be determined in the same way 
as a one-step bank for the given values of S m i n , Qd.min and Qo.max- The function of a two-step hierarchy is, therefore, 
limited to providing an estimate of the location of the global maximum. For this it utilizes the information that the 
occurrence of a (high) threshold crossing must generate in templates that are relatively far away from it. 

We give here a brief review of the algorithm used to set up a two-step hierarchical search in the Newtonian case. 
Let the spacing between consecutive templates (in terms of To) in the one-step template bank, for given S'min, <3d,min 

and Qo,max, be S 2 . Then the first stage spacing Si is taken to be Si — k5%, k = 2, 3, For each Si, the first stage 

threshold r/i was kept such that a signal with strength S'min, lying in the middle of two consecutive templates, has 
a detection probability Qd,min- The average number of false crossings among the first stage templates can then be 
computed which, given the number of second stage templates to employ around each first stage crossing, in turn 
allows the overall average computation cost to be calculated. The minimum of this cost was then found as a function 
of k. 

In the present case, a similar approach can be followed to space the first stage templates as was done for the 
Newtonian case. Here, there will be two spacings to fix, namely, along the minor and major axes of the one-step unit 
cells, and these can be chosen as integral multiples of the corresponding one step spacings l\ and l 2 . For convenience 
in the following, we denote a first stage unit cell at the location 9 a as Uj{ki, k 2 , 9 a ), where the sides of the unit cells 
have lengths ki x li and k 2 x l 2 . There are, however, a few complications that arise in this approach. First, if a first 
stage unit cell Ui(ki,k 2 , 9 a ) satisfies CI, it is not implied that a unit cell Ui(ki,k 2 , 9b) at a different location will 
also do so. This is because the dimensions of a first stage unit cell would be quite large and Figs. |l3|, |IJ show that 
the error in detection probability rises with an increase in the dimensions. Thus the same (k±, k 2 ) at two different 
locations would lead to different values of 771 since the first stage threshold is determined by detection probability. 
We take this effect into account as follows. For each (k±, k 2 ), we take two widely separated locations 8 a and Ob and 
compute the thresholds r]\ a and r)ib that are required to make both Ui(ki,k 2 , 9 a ) and Ui(ki,k 2 , 9b) satisfy CI. We 
then choose the minimum among these two as the first stage threshold 771 . 

The second complication is the boundary near vertex A of the space of interest which also disallows the same 
(fci, k 2 ), as in the broader parts of the space of interest, from being used in this region. Recall that the one step 
template grid in this region was constructed out of single unit cells "stacked" vertically (see Fig. |l5|). Hence, if k 2 > I, 
extra second stage templates would be required in the region to the right of AB which would increase the number of 
one-step templates without adding significantly to the range of binary masses being detected. However note that if 
Ui(ki,k 2 , 9 a ), for k 2 > 1, satisfies CI for some threshold 771, then so would Ui(k±, 1, 9 a ), since the templates at the 
vertices will now be closer to all the signals in the cell's interior. Therefore, while calculating the number of first stage 
templates in this region, we simply divide the number of one step templates by k\ . 

Finally, the number of second stage templates that would be employed per first stage crossing will now depend on 
whether the crossing occurs in the narrow region near vertex A or in the broader part of the space of interest. Let 
this number be n. In the broader part of the space of interest, 

n = 4 (h - l)(k 2 - 1) + 2 (ki - 1) + 2 (k 2 - 1) . (93) 

Now, the minimum of the computational cost of a two-step search occurs when the number of false crossings in the 
first stage becomes ~ I. But most of the first stage templates will be located in the broader part of the space of 
interest and, hence, most of the false crossings will also occur in this region. This implies that n will be as given in 



Eq. (93) for most cases. We, therefore, take n to be the above for all crossings. Note that, this assumption would 
lead to an overestimate of the computational requirements for the first stage and is, thus, "safe" in this sense. 

Let the number of false crossings for a given input data segment be n c . Then the total number of templates which 
will be employed for that data segment would be + n c n, where is the total number of first stage templates. 
In the presence of a signal, there would be an extra term of n in the above sum but since the event rate of signals 
is expected to be quite low, this term can be neglected. If we assume that the first stage rectified output are all 
statistically independent of each other, as would be the case if they are spaced widely apart, then the average number 
of false crossings nf would be, 

nf = n c x Qo(m) ■ (94) 
Qoivi) is the probability of at least one crossing of ?/i in a single rectified output (see Eq. (|63|)), 
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Qo(x) = 1 - exp[-T> s exp(-a; 2 /2)] , (95) 



where v s is the sampling rate and Tp was defined in Section II B to be the padding in the time series of the template 
wave form having the longest duration. Note that since rji would be ~ 6.0, the effective statistical independence 
of rectified output samples may be less and e should be less than unity. However, keeping e — 1 leads to an over 
estimation of n* v and, hence, an under estimation of the computational advantage of a two step search. 
The average total computational cost for a two-step search would be, 

= n ( p + ra* v x n . (96) 
In order to compare the performance of a two-step search with the corresponding one-step search (that is, for the same 

*Jminj Qd, min 

and Qo.max), we use the computational powers required for implementing the two strategies on-line (that 
is, the input data should be processed in the same time as required in its collection). The number of floating point 
operations required in a one-step search to process T sec of data can be estimated as follows: (i) The number of flop 
involved in the Discrete Fourier transform (DFT) x of the detector output time series x would be 3 N log 2 N where 
N = v s T. However, this transform needs to be computed only once and, thus, does not contribute significantly to 
the total computational cost, (ii) For each template location 6 a , two correlations would be required, namely, with the 
quadrature components q and 5^/2- This involves computing the product of x with the DFTs of q and q^/ 2 followed 
by an inverse DFT for each of the resulting series. Hence, 2 N + 6 N log 2 N flop will be required here, (iii) Each of 
these transformed series would then have to be squared and added but only the first Tp sec of each series is required. 
This, thus, leads to 3Tp v s flop. Thus, the total number of operations, -/V nop involved in a one step search is, 

jV flop = jyW x (2 N + 6 N log 2 N + 3 T P v s ) , (97) 
For an on-line one-step search, iVfl op operations would have to be performed in Tp sec. Thus, 

C£L = >< 10~ 9 Gnops , (98) 

A similar estimate for an on-line two step search leads to an average computational requirement of, 

N {2) 

C£L = x ID" 9 Gnops (99) 

where, 

N^ p = X (2 JV + 6 N log 2 N + 3 T P v s ) . (100) 

We call the quantity Conline/^onliae as t ne computational advantage C ga i n of a two-step search. This is the factor by 
which a two-step hierarchical search would be faster than the corresponding one-step search in an on-line detection. 

We now present our results in the form of Table |[ for initial LIGO, and Table || for advanced LIGO. The value of 
5 m i n for each table has been taken sufficiently low so that the resulting one-step unit cells obtained are small. It was 
shown in MD96 that, for a given number of templates, the one-step threshold is almost independent of T for low false 
alarms. This implies that the unit cell dimensions will also be independent of T (the variation of the threshold in the 
advanced LIGO case is larger but it is still negligible). Therefore, the values of It and Z 2 , for the one-step unit cell, 
are given in the caption of each table. These values are for a unit cell located at (1.3, 25.0) sec for the case of initial 
LIGO and (13.0, 1000.0) sec for the advanced LIGO. The values of the one-step threshold (which is the second stage 
threshold rf^ for the two-step search) and the total number of one-step templates (obtained by taking the variation 
of unit cell areas into account) is also given. 

The first column in each table is the value of T. Since the sampling rate used in our calculation is 2 11 = 2048 Hz 
and an FFT is most efficient when the number of samples is a power of two, we choose T to be a power of 2 also. The 
second and third columns are the values of fci and fc 2 at which the average computational cost of the two-step search 
is minimised. The fourth column is the corresponding first stage threshold n^ and the fifth and sixth columns are 
the corresponding values of n^ v and nf. We have kept only the integral part of n^ v and nf and, therefore, n^ v = 
means that n^ v ~ 1 or less. The seventh column is the computational power required for an on-line two-step search 
followed by the computational power required for an on-line one-step search in the eigth column. The last column 

lists Cgain- 
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Even though we have used large values of T, especially in Table |ll| such values would be difficult to use in a practical 
implementation because of memory restrictions. We have used these values only to show the existence of a minimum 
in the computational power requirement as a function of T p|| . It should be noted here that for the case of advanced 
LIGO, the storage of the pre-computed template wave forms is also a significant problem. For instance, even if we 
consider the average duration of templates in the advanced LIGO case to be ~ 100 sec, the amount of storage required 
for all the ~6x 10 5 quadrature Fourier transforms would be ~ 100 x 2048 x 6 x x8/10 = 983 GB (assuming that 
each sample value requires 8 Bytes of storage). This is a very low bound since the duration of a significant number 
of templates will be much larger. 

The results obtained above can be checked approximately as follows. The one-step template placement criterion 
of 0] requires the templates to be placed such that, for any signal, H = 0.97 in at least one nearby template. Then 
the number of one-step templates N^, would be the area A of the space of interest divided by the area of the 0.97 
contour. For advanced LIGO, Nj, — 20389.5/0.04 ~ 509739. Thus, the threshold rj 2 required, for a false alarm rate 
of 1 false event/y, would be rf 2 > = 8.722 for T = 8192.0 sec. For the detection probability formula used in MD96, 
it was found that the minimum observed strength required for a signal so that its detection probability be 0.95 is 
•Sobs ~ V + 0.67. The actual strength should, therefore, be S m i n = 5 o bs/0.97 = 9.682. Roughly speaking, the 
decrease in nf , with an increase in k\ and k% , is halted when becomes of order unity. Assuming that the number 
of first stage templates that is finally obtained is ~ 10 4 , it would imply that, for the above value of T, r^ 1 ) » 7.026. 
For a detection probability of 0.95 in the first stage, therefore, the value of S bs = S m i n H' should be 7.697, where TL' is 
the value of the intrinsic ambiguity in the middle of the sides of a first stage unit cell, i.e., HI = H(0 a , 8 a + kiliei a /2). 
The quantity ki can then be calculated as the ratio of the dimension of the TL' contour along ei a to (0.97, 9 a ). From 
the above, Hi = 0.79 which gives k\ — 7.67, ki — 4.35 (we have allowed ki to be non-integral here). These values are 
about the same as those in Tables [0. However, this approximation is crude in many ways and can only serve as an 
indicator for the kind of values one may get for ki. 

The savings in computational requirements achieved by a two-step search can be more than what is obtained here 
if the first stage template grid is rotated relative to the second stage grid. This is because of the shear of the contours. 
In the argument given above, the quantities k\ and ki were obtained as the ratios of Zi(0.97, a ) and the corresponding 
dimension of the lower level contour TL 1 . However, the direction in which the TL' contour is most elongated is different 
from that of the eigenvectors ei . If the first stage grid were oriented along the direction of maximal elongation of TL' 1 
the first stage unit cell may turn out to be larger. However, the calculation of the number of first stage templates as 
well as the number, n, of second stage templates would be more involved in such a case. We postpone an investigation 
of this problem to a later work. 

VI. CONCLUSIONS 

We have investigated the performance of a two step hierarchical search for the detection of gravitational wave 
signals emitted during the inspiral of a compact binary. This work extends the investigations of MD96 || to the more 
realistic case of zero spin post 15 -Newtonian template and signal wave forms. 

As in MD96, we find that a two-step search brings about a significant reduction in computational requirements. 
For the case of (i) initial LIGO noise p.s.d., a two-step search is ~ 27.0 times faster than the corresponding one-step 
search and (ii) for the advanced LIGO noise p.s.d., a two-step search is ~ 23.0 times faster than the corresponding 
one-step search. The range used for the masses mi and m-i is 0.5 < m\ < 30.0 M©, 0.5 < m-i < 30.0 M©. 

In the analysis of MD96, the dominant problem was the calculation of detection probability in the presence of 
strong statistical correlations between rectified output samples. A solution to this problem was found in this paper 
in the form of a semi-analytic method that reproduces the exact Monte Carlo estimates quite well. It is also shown 
here that statistical correlations are unimportant for the calculation of false alarm probability when the threshold is 
kept sufficiently high. Therefore, the effective sampling rate used in MD96 is not required. 

Though the issues of detection and false alarm probabilities have been addressed satisfactorily here, some new 
problems crop up in the present analysis. Namely, the (i) location dependence of the intrinsic ambiguity function and 
(ii) the non-trivial shape of the boundary of the space of interest. Both these problems were dealt with by making 
some approximations. The location dependence of the intrinsic ambiguity function seems weak enough, at least in 
the case of advanced LIGO, for us to assume that the grid of one-step templates will be an "adiabatically" changing 
regular grid. This allows us to approximately take the effect of variations in the area of the contours into account. 
The non-trivial boundary has a significant effect only near one of the vertices (vertex A of Fig. |l|). We take this effect 
into account by placing a single "stack" of unit cells in this region. 
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The results of this paper show that the use of hierarchical methods of detection can be very useful for the case of 
coalescing binary signals and provide a strong motivation for more detailed investigations. Such methods would be 
indispensable if the number of signal parameters required becomes large. For instance, if the orbital and total angular 
momenta of the binary are misaligned, there would be significant modulations of the phase and amplitude which can 
reduce the signal to noise ratio if these effects are neglected in the template family. For such signals, a template family 
with a larger number of parameters may be required. 

Many other hierarchical strategies are also conceivable and it remains to be seen whether they can be more effective 
than the two-step search analysed here. For instance, one obvious strategy is to use a lower order template family as 
the first stage of the search and use the true wave forms, having a larger number of parameters, as the second stage. 
It is not enough, though, to only provide estimates of their performance since at some stage such strategies need to 
be implemented in practice and, as seen in this paper, the details of the implementation can also be an involved issue. 
Also, the robustness of the placement configuration against changes in the noise power spectral density needs to be 
investigated. The efficacy of hierarchical methods (not necessarily a two-step search) should also be investigated for 
the detection of continuous wave sources where the estimated computational requirements are extremely large and 
far beyond presently available computing power. Further investigations in this direction are in progress. 
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APPENDIX A: THE BIVARIATE PROBABILITY DENSITY P Zl ,z 2 AND Z\ Z 2 

Here, we outline the steps in the derivation of Eq. (p5[). The algebraic manipulations were performed using 
mat hematic A . First, the general expression for the joint bivariate probability density is derived without assum- 
ing the mean values of the Gaussian components to be zero. Let the bivariate cumulative distribution function of 
Z x = [Xl + Xl] 1 / 2 and Z 2 = [Yf+Yg] 1 / 2 be F Zl ,z 2 ( z i, z i), where (X 1 , X 2 ,Y 1 ,Y 2 ) is a set of jointly Gaussian random 
variables with a covariance matrix give in Eq. ( |48| ) and mean values X x = fj,i, X 2 = /j,2, Y\ = ui, Y% = v^. Changing 
the variables of integration to X\ = Rcostfi, X 2 = Rsm.<f), Y% — Qcos-0 and Y 2 — Qsixitp, we get, 
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Eq. (Al) can be rewritten as, 
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where, 
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E = [(v\ - r/ii + s/i 2 ) 2 + (v 2 - r^ 2 - s^if] 1/2 

1 /2 

D = [(/i2 - rv 2 + svx) 2 + {fix - rv Y - sv 2 ) 2 ] 
Xi = arctan 
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The integral over (f> can be performed to yield, 
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where, Iq(x) is the modified Bessel function of the first kind of order zero. The probability density function, Pz 1 ,z 2 > 
can be obtained now as, 
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In the absence of a signal, /ii = /X2 = v\ = v 2 = and the joint probability density reduces to, 
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Thus, the correlation TTv can be obtained as, 
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The above double integral is solved in (24j from which we get 
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where E is the complete elliptic integral of the second kind and K is the complete elliptic integral of the first kind. 
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TABLE I. 
h = 0.022 sec 


Minimum C™_ 
, h = 0.144 sec, 


■line 
AT*, = 


as a function of T for : 
= 13279. 


Smin — 


9.98, Z max = 140.482 sec, 


Qd,m*n = 0.95, 7? (2) = 


8.314, 


T(sec) 


fci k-i 






< 


< 


C%LJP flops) 


C^jq flops) 


Cgain 


256.0 

512.0 

1024.0 

2048.0 

4096.0 

8192.0 


8 9 
8 6 
8 5 
8 4 
8 4 
8 3 




6.056 
6.283 
6.484 
6.649 
6.649 
6.866 







1 




360 
441 
490 
620 
682 
733 


0.192 
0.155 
0.152 
0.187 
0.207 
0.228 


7.07 
4.65 
4.12 
3.99 
4.02 
4.12 


36.82 
30.00 
27.11 
21.34 
19.42 
18.07 


TABLE II. 
h = 0.116 sec 


Minimum C™_ line 
, h = 0.560 sec, iV£ = 


as a function of T for : 
= 300796. 


Smin 


10.34, £, max = 5621.51 sec, 


Qd, mln = 0.95, ?7 (2) = 


8.658, 


T(sec) 


fci &2 






nf 


nT 


C^JG flops) 


d£LjG flops) 


Cgain 


8192.0 

16384.0 

32768.0 

65536.0 

131072.0 


5 9 
4 7 
4 7 
4 7 
4 6 




6.649 
7.002 
7.002 
7.002 
7.060 


11 

6 
16 
35 
56 


10188 
13490 
14390 
16100 
18935 


9.771 
6.476 
5.709 
6.014 
7.004 


288.48 
144.39 
119.34 
112.36 
111.27 


29.52 
22.30 
20.90 
18.68 
15.89 
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FIG. 1. The space of interest for the case of initial LIGO noise power spectral density (f a = 40.0 Hz). The vertices of the 
space of interest correspond to the following points in the (mi, 7712) plane : A corresponds to (0.5, 0.5) M©, B to (30.0, 0.5) M© 
and C to (30.0, 30.0) M©. 

FIG. 2. The space of interest for the case of advanced LIGO noise power spectral density (f a = 10.0 Hz). The vertices of the 
space of interest correspond to the following points in the (mi, m 2 ) plane : A corresponds to (0.5, 0.5) M©, B to (30.0, 0.5) M© 
and C to (30.0, 30.0) M©. 

FIG. 3. The relative error in detection probability as obtained using the multivariate Gaussian approximation and 
as obtained by performing an exact simulation. Each figure shows the relative error for three values of signal strength, 
S — 8.0 (solid), 9.0 (dotted), 10.0 (dashed). As expected, the error decreases for larger signal strengths. The top left figure 
shows the locations of the templates (crosses), used in the calculation of the detection probabilities, and the signal loca- 
tions (filled circles). The basic unit cell which is composed of templates # 1, 2, 3 and 4, is oriented along the eigenvectors of 
the Hessian at the location of template # 1. Templates # 5 and # 6 are included in the calculation to take into account any 
possible contribution that they may provide, because of the shear of contours, to the detection probabilities of signal # 3 and 
# 2 respectively. 

FIG. 4. The contours of the intrinsic ambiguity function TL{9 a , 9b) for initial LIGO. In this figure, 9 a is kept fixed at 
8 a = (1-3, 25.0) sec and 9b is varied. Also shown are the semi-minor and semi-major axes of the 0.97 contour as calculated 
from the Hessian H(# a ). The axes do not appear at a right angle to one another because the axes scales are different. 

FIG. 5. The upper figure shows the contours of ^(0.97, 9 a ) and the lower figure shows the contours of Zi (0.97, 9 a ), on the 
space of interest for initial LIGO. 

FIG. 6. The upper figure shows the area, 7rZi(0.97, 9 a ) Z2(0.97, 9 a ) of the 0.97 contour of the intrinsic ambiguity function H 
and the lower figure shows the contours of an, the angle (in degrees) between ei a and the T1.5 axis, on the space of interest for 
initial LIGO. 

FIG. 7. The upper figure shows the contours of ^(0.97, 9 a ) and the lower figure shows the contours of Zi(0.97, 9 a ), on the 
space of interest for advanced LIGO. 

FIG. 8. The upper figure shows the area, 7rZi(0.97, 9 a ) ^(0.97, 9 a ) of the 0.97 contour of the intrinsic ambiguity function TL 
and the lower figure shows contours of on, the angle (in degrees) between ei a and the T1.5 axis, on the space of interest for 
advanced LIGO. 

FIG. 9. The parameters of a unit cell which is a parallelogram. Whenever the sides of the unit cell are assumed to be along 
the eigenvectors of the Hessian, l\ is taken to be along the semi-minor axis while I2 is taken along the semi-major axis. 

FIG. 10. In the uppermost figure we show the detection probability of signals in the interior of a rectangular unit cell which 
is oriented along the T1.5 (along the x-axis in the figure) and to (along the y-axis) axes with the top left vertex at (1.1, 25.0) sec. 
The length h of the side along the to axis is 0.05 sec while the length I2 of the other side is 0.150 sec. The threshold and signal 
strengths were chosen (arbitrarily) as 77 = 8.0 and S = 9.0. In the middle figure the contours of this detection probability map 
are superimposed on some of the contours (dashed) of H(9 a , 9b) with 9 a = (1.1, 25.0) sec. In the lowermost figure, we show 
the detection probability map for the same unit cell but now with l\ and I2 oriented along the eigenvectors ei a and eia- The 
threshold and signal strength are the same as in the figures above it. The values along the x-axis and y-axis are the serial 
numbers of the grid points. 

FIG. 11. The quantities Si, 82 and S a (see Eq. (83) and Eq. (84)) for the case of initial LIGO. Each row of figures corresponds 
to a fixed value of to. For (i) the first row, to = 60.0 sec (ii) the second row, to = 30.0 sec and (iii) for the third row, to = 10.0 sec. 
The first column corresponds to 5i, the second corresponds to £2 and the third to 5 a . The x-axis is the T1.5 axis. 
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FIG. 12. The quantities 61, S 2 and 5 a (see Eq. (64) and Eq. (65)) for the case of advanced LIGO. Each row of figures 
corresponds to a fixed value of To. For (i) the first row, To = 3000.0 sec (ii) the second row, tq — 2000.0 sec and (iii) for the 
third row, To = 800.0 sec. The first column corresponds to 5i, the second corresponds to 62 and the third to S a . The abscissa 
is the T1.5 axis. 

FIG. 13. The relative error in detection probability for corresponding signals in two widely separated unit cells. Let the 
locations of the unit cells be 6 a and 8b- In this figure, we consider the case of initial LIGO and place the top left templates of 
the unit cells at 8 a — (1.3, 50.0) sec and 9b = (1.5, 10.0) sec. The value used for the signal strength S is shown at the top of 
each plot. The effect of plunge cutoff has been incorporated in the calculations. 

FIG. 14. The relative error in detection probability for corresponding signals in two widely separated unit cells. Let the 
locations of the unit cells be 8 a and 9b- In this figure, we consider the case of advanced LIGO and place the top left templates 
of the unit cells at 6 a = (13.0, 2000.0) sec and 8 b = (15.0, 400.0) sec. The value used for the signal strength S is shown at the 
top of each plot. The effect of plunge cutoff has been incorporated in the calculations. 

FIG. 15. A schematic illustration of a quasi-regular grid of unit cells near the vertex A of the space of interest (initial LIGO). 
The lengths used for the sides of the unit cells are h — 0.02 sec and fo = 0.120 sec. These lengths have been chosen arbitrarily 
but represent typical values obtained in a one-step search. The boundary of the space of interest is shown by the lighter lines. 
The x-axis represents T1.5 and to lies along the y-axisn (both units are in seconds). The top left corner of each unit cell is 
placed on the left most boundary which is the image of the principal diagonal in the (mi, m^) plane. 

FIG. 16. The number of templates for the case of initial LIGO as a function of the unit cell dimensions h and fo. The solid 
contours are obtained by using the algorithm that takes the variation of unit cell areas into account. The dashed contours are 
for the values obtained by simply dividing the area of the space of interest by h x fo. 

FIG. 17. The number of templates for the case of advanced LIGO as a function of the unit cell dimensions h and fo. The 
solid contours are obtained by using the algorithm that takes the variation of unit cell areas into account. The dashed contours 
are for the values obtained by simply dividing the area of the space of interest by h x fo. 
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